!----------------------------------------------------------------------!
! Optional advection scheme for CCATT-BRAMS/BRAMS models version 4.2+  ! 
! Based on Walcek, 2000 (JGR) and Walcek and Aleksic, 1998 (ATENV).    ! 
! The scheme is highly conservative, monotonic and keeps mass mixing   ! 
! ratio positive definite. 					       ! 
! Implemented by Saulo Freitas (saulo.freitas@cptec.inpe.br) @ Jun/2009! 
!----------------------------------------------------------------------!

MODULE monotonic_adv
  
  USE node_mod, only:        &
          	      ibcon, &  !INTENT(IN)
          	      mynum, &  !INTENT(IN)
          	      nodei0,   &  !INTENT(IN)
          	      nodej0,   &  !INTENT(IN)
          	      nodemyp,  &  !INTENT(IN)
          	      nodemxp,  &  !INTENT(IN)
          	      nodemzp   	!INTENT(IN)
  USE mem_grid, ONLY:        &
  	    	     dtlt,   & !INTENT(IN)\
		     time,   &
  	    	     ngrids, & !INTENT(IN)
  	    	     ngrid,  & !INTENT(IN)
	    	     dzt,    & !INTENT(IN)
	    	     dztn,   & !INTENT(IN)
            	     grid_g, & !INTENT(IN)	 
	    	     grid_g, & !INTENT(IN) 
	    	     naddsc, & !INTENT(IN) 
	             hw4       !INTENT(IN)	   
  
  USE mem_basic, ONLY: basic_g  !INTENT(IN) 
  
  USE micphys    ,  ONLY: level !INTENT(IN) 
  
  USE rconstants ,  ONLY: cp,p00,cv,rgas,cpi   !INTENT(IN)

  !lfr use mem_aer1, only : aerosol,num_scalar_aer_1st   !INTENT(IN)
  use mem_chem1, only : nspecies_transported        !INTENT(IN)

  use mem_scratch, only  : scratch  ! only scr1, inout
  
  use var_tables, only : scalar_tab & ! (var_p = IN, var_t = INOUT) 
                        ,num_scalar   ! (IN) 
			
  


  !LFR -> To test
  !USE stopwatch


  IMPLICIT NONE
  include "mpif.h"
  
  INTEGER:: advmnt

  INTEGER , PARAMETER :: ON=1,OFF=0
  
  INTEGER, PARAMETER :: GhostZoneLength=15

  INTEGER , PARAMETER :: use_true_density  = 0 ! 0= OFF, 1=ON
  
  !- for theoretical experiments
  INTEGER , PARAMETER :: theor_wind = 0        ! 0= OFF, 1=ON

  INTEGER           :: mnt_adv_jnitialized=0
  
  real, parameter :: c1 = cv/rgas, c2 = p00/rgas !c2 = p00*(cpi**c1)/rgas

  TYPE advmnt_vars 
     REAL,POINTER,DIMENSION  (:,:,:)  :: u3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: v3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: w3d 

     REAL,POINTER,DIMENSION  (:,:,:)  :: vc3d_jn 
     REAL,POINTER,DIMENSION  (:,:,:)  :: vc3d_out 
     REAL,POINTER,DIMENSION  (:,:,:)  :: vc3d_x 
     REAL,POINTER,DIMENSION  (:,:,:)  :: vc3d_y 

     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3du 
     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3dv 
     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3dw 

     REAL,POINTER,DIMENSION  (:,:,:)  :: den0_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: den1_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: den2_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: den3_3d 

     REAL,POINTER,DIMENSION  (:,:)    :: dxtW 
     REAL,POINTER,DIMENSION  (:,:)    :: dytW 
     REAL,POINTER,DIMENSION  (:)      :: dztW 

  END TYPE advmnt_vars

  TYPE(advmnt_vars), ALLOCATABLE,DIMENSION(:)   :: advmnt_g,advmntNew_g
    
  type communicate
     integer,pointer,DIMENSION(:) :: tag
     integer,pointer,DIMENSION(:) :: length
     integer,pointer,DIMENSION(:) :: procSend
     integer,pointer,DIMENSION(:) :: iaSend
     integer,pointer,DIMENSION(:) :: izSend
     integer,pointer,DIMENSION(:) :: jaSend
     integer,pointer,DIMENSION(:) :: jzSend
     integer,pointer,DIMENSION(:) :: procRecv
     integer,pointer,DIMENSION(:) :: iaRecv
     integer,pointer,DIMENSION(:) :: izRecv
     integer,pointer,DIMENSION(:) :: jaRecv
     integer,pointer,DIMENSION(:) :: jzRecv
  end type communicate
  TYPE(communicate),allocatable,DIMENSION(:) :: comm_i,comm_j

  PUBLIC :: advmnt_driver  ! Subroutine
  PUBLIC :: procMap

  INTEGER, PARAMETER :: nvariables=2

  INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: procMap
  INTEGER, ALLOCATABLE, DIMENSION(:,:,:,:) :: itag1,itag2
  INTEGER, ALLOCATABLE, DIMENSION(:) :: handle

  LOGICAL,DIMENSION(:),ALLOCATABLE :: west,east,north,south
  INTEGER,ALLOCATABLE,DIMENSION(:,:,:,:,:) :: lineTag
  INTEGER,ALLOCATABLE,DIMENSION(:,:,:,:,:) :: columnTag

  !nvariables=
  INTEGER, PARAMETER :: v_flux=1
  INTEGER, PARAMETER :: v_jmxmn=2
  
  
  INTEGER :: numberOfProcessors

  INTEGER,ALLOCATABLE :: procRecv_i(:),      procRecv_xi(:)
  INTEGER,ALLOCATABLE :: tagRecv_i(:),	     tagRecv_xi(:)
  INTEGER,ALLOCATABLE :: iaRecv_i(:),	     iaRecv_xi(:)
  INTEGER,ALLOCATABLE :: izRecv_i(:),	     izRecv_xi(:)
  INTEGER,ALLOCATABLE :: jaRecv_i(:),	     jaRecv_xi(:)
  INTEGER,ALLOCATABLE :: jzRecv_i(:)	,    jzRecv_xi(:)    	   
  INTEGER,ALLOCATABLE :: bufRecvLength_i(:), bufRecvLength_xi(:)
  INTEGER,ALLOCATABLE :: bufRecvStart_i(:),  bufRecvStart_xi(:)
  INTEGER,ALLOCATABLE :: procSend_i(:),      procSend_xi(:)
  INTEGER,ALLOCATABLE :: tagSend_i(:),	     tagSend_xi(:)
  INTEGER,ALLOCATABLE :: iaSend_i(:),	     iaSend_xi(:)
  INTEGER,ALLOCATABLE :: izSend_i(:),	     izSend_xi(:)
  INTEGER,ALLOCATABLE :: jaSend_i(:),	     jaSend_xi(:)
  INTEGER,ALLOCATABLE :: jzSend_i(:),	     jzSend_xi(:)
  INTEGER,ALLOCATABLE :: bufSendLength_i(:), bufSendLength_xi(:)
  INTEGER,ALLOCATABLE :: bufSendStart_i(:),  bufSendStart_xi(:)
  INTEGER,ALLOCATABLE :: procRecv_j(:),      procRecv_xj(:)
  INTEGER,ALLOCATABLE :: tagRecv_j(:),	     tagRecv_xj(:)
  INTEGER,ALLOCATABLE :: iaRecv_j(:),	     iaRecv_xj(:)
  INTEGER,ALLOCATABLE :: izRecv_j(:),	     izRecv_xj(:)
  INTEGER,ALLOCATABLE :: jaRecv_j(:),	     jaRecv_xj(:)
  INTEGER,ALLOCATABLE :: jzRecv_j(:)	,    jzRecv_xj(:)    
  INTEGER,ALLOCATABLE :: bufRecvLength_j(:), bufRecvLength_xj(:)
  INTEGER,ALLOCATABLE :: bufRecvStart_j(:),  bufRecvStart_xj(:)
  INTEGER,ALLOCATABLE :: procSend_j(:),      procSend_xj(:)
  INTEGER,ALLOCATABLE :: tagSend_j(:),	     tagSend_xj(:)
  INTEGER,ALLOCATABLE :: iaSend_j(:),	     iaSend_xj(:)
  INTEGER,ALLOCATABLE :: izSend_j(:),	     izSend_xj(:)
  INTEGER,ALLOCATABLE :: jaSend_j(:),	     jaSend_xj(:)
  INTEGER,ALLOCATABLE :: jzSend_j(:),	     jzSend_xj(:)
  INTEGER,ALLOCATABLE :: bufSendLength_j(:), bufSendLength_xj(:)
  INTEGER,ALLOCATABLE :: bufSendStart_j(:),  bufSendStart_xj(:)


  INTERFACE get
     MODULE procedure &
          get_real, &
          get_int, &
          get_logical
  END INTERFACE

  INTERFACE send
     MODULE procedure &
          send_real, &
          send_int, &
          send_logical
  END INTERFACE

  !type (watchtype) :: w(4)   ! This declares w to be a watch
  !type (watchtype) :: w1(6)   ! This declares w to be a watch
  
  INTEGER :: nSend_i,nSend_j,nRecv_i,nRecv_j
  INTEGER,PARAMETER :: bigdump=1
  REAL,ALLOCATABLE,DIMENSION(:,:,:) :: buffcomm
  INTEGER,ALLOCATABLE,DIMENSION(:) :: newM2, newM3, newIa, newJa
  INTEGER :: nRec_i,nSnd_i,nRec_j,nSnd_j
  INTEGER :: bufSendTotalLength_i,bufSendTotalLength_j
  INTEGER :: bufREcvTotalLength_i,bufRecvTotalLength_j

CONTAINS
  !----------------------------------------------------
             !call ADVECTc('V' ,mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum)
  SUBROUTINE advmnt_driver(varn,m1 ,m2 ,m3 ,ia,iz,ja,jz,izu,jzv,mynum)
  use node_mod, only : nodeibcon,nodeia,nodeiz,nodeja,nodejz,ipaths,iget_paths
    IMPLICIT NONE
    INTEGER , INTENT(IN) :: m1
    INTEGER , INTENT(IN) :: m2
    INTEGER , INTENT(IN) :: m3
    INTEGER , INTENT(IN) :: ia
    INTEGER , INTENT(IN) :: iz
    INTEGER , INTENT(IN) :: ja
    INTEGER , INTENT(IN) :: jz
    INTEGER , INTENT(IN) :: izu
    INTEGER , INTENT(IN) :: jzv
    INTEGER , INTENT(IN) :: mynum
    character(len=*),intent(IN) :: varn
    
    !--- local vars
    integer n,ng,mxyzp,i,j,procfile
    real, pointer :: scalarp, scalart
    
    character(4) :: ctime
        
    INTEGER :: iteste !lfr
    !PRINT *,'No driver da advect ',mynum;CALL flush(6)

    !call create_watch(w,name=(/"Waiting MPI "/))    ! Watches must be created before they are used
    !-scratch arrays initialization  and air density for using
    !- basic state densities
    
    !PRINT *, 'In advection........time=',time;CALL flush(6)
    
    iteste=0 !lfr
    
    IF(mnt_adv_jnitialized == OFF) THEN
       PRINT *,'Allocating Communication info variables in proc. ',mynum; CALL flush(6)
       CALL CommAlloc(mynum,ngrid,m1,m2,m3,ia,iz,ja,jz)

       PRINT *,'Fill local Communication info variables in proc. ',mynum; CALL flush(6)
       CALL BuildComm(mynum,ngrid,m1,nRec_i,nSnd_i,nRec_j,nSnd_j, &
                        bufSendTotalLength_i,bufSendTotalLength_j, &
			bufRecvTotalLength_i,bufRecvTotalLength_j)

       IF(bigdump==1) CALL dumpBorders('Borders-','i',m1, m2, m3,mynum, &
	   nRec_i, procRecv_i, tagRecv_i, &
	   iaRecv_i, izRecv_i, jaRecv_i, jzRecv_i, &
	   bufRecvStart_i, bufRecvLength_i, bufRecvTotalLength_i, &
	   nSnd_i, procSend_i, tagSend_i, &
	   iaSend_i, izSend_i, jaSend_i, jzSend_i, &
	   bufSendStart_i, bufSendLength_i, bufSendTotalLength_i)     

       IF(bigdump==1) CALL dumpBorders('Borders-','j',m1, m2, m3,mynum, &
	   nRec_j, procRecv_j, tagRecv_j, &
	   iaRecv_j, izRecv_j, jaRecv_j, jzRecv_j, &
	   bufRecvStart_j, bufRecvLength_j, bufRecvTotalLength_j, &
	   nSnd_j, procSend_j, tagSend_j, &
	   iaSend_j, izSend_j, jaSend_j, jzSend_j, &
	   bufSendStart_j, bufSendLength_j, bufSendTotalLength_j)     


       CALL initialize_advmnt(ngrids,nodemzp(mynum,:), &
                              nodemxp(mynum,:),nodemyp(mynum,:))


       do ng=1,ngrids
         CALL initialize_grid_spacings(ng,nodemzp(mynum,ng), &
	               nodemxp(mynum,ng),nodemyp(mynum,ng) &
                      ,grid_g(ng)%dxt	    &
                      ,grid_g(ng)%dyt	    &
 		      ,grid_g(ng)%fmapt     &
		      ,grid_g(ng)%rtgt      &
!
    		      ,advmnt_g(ng)%dxtW &
    		      ,advmnt_g(ng)%dytW &
		      ,advmnt_g(ng)%dztW )
       enddo
       
       if(use_true_density == OFF) then 
          do ng=1,ngrids
	    CALL initialize_densities(nodemzp(mynum,ng),&
	               nodemxp(mynum,ng),nodemyp(mynum,ng) &
                      , basic_g(ng)%dn0     &
		      , basic_g(ng)%dn0u    &
                      , basic_g(ng)%dn0v    &
                      ,advmnt_g(ng)%dd0_3d  &
    		      ,advmnt_g(ng)%dd0_3du &
    		      ,advmnt_g(ng)%dd0_3dv &
		      ,advmnt_g(ng)%dd0_3dw )
	   enddo
       endif
       
       mnt_adv_jnitialized= ON    
    END IF

    mxyzp=m1*m2*m3
    
    !- Advect  U, V, and W
    if (varn .eq. 'V' .or. varn .eq. 'ALL') then
      stop 'not using mnt to advect u,v,w'
    endif

    !- Advect  scalars
    
    !- get actual air densities, if using them instead of basic state fields
    if(use_true_density == ON) then
     call get_true_densities(m1,m2,m3,level &
                      , basic_g(ngrid)%rtp     &
		      , basic_g(ngrid)%rv      &
                      , basic_g(ngrid)%pp      &
		      , basic_g(ngrid)%pi0     &
                      , basic_g(ngrid)%theta   &
                      ,advmnt_g(ngrid)%dd0_3d  &
    		      ,advmnt_g(ngrid)%dd0_3du &
    		      ,advmnt_g(ngrid)%dd0_3dv &
		      ,advmnt_g(ngrid)%dd0_3dw )
    
    endif

 
    !- prepare wind velocities including map factors
    call prepare_winds(dtlt,m1,m2,m3,ia,iz,ja,jz     &
                      ,basic_g(ngrid)%uc  &
		      ,basic_g(ngrid)%up  &
                      ,basic_g(ngrid)%vc  &
                      ,basic_g(ngrid)%vp  &
		      ,basic_g(ngrid)%wc  &
                      ,basic_g(ngrid)%wp  &
!        
 		      ,grid_g(ngrid)%fmapui &
		      ,grid_g(ngrid)%fmapvi &
		      ,grid_g(ngrid)%rtgt   &
                      ,grid_g(ngrid)%rtgu   &
		      ,grid_g(ngrid)%rtgv   &
		      ,grid_g(ngrid)%f13t   &
                      ,grid_g(ngrid)%f23t   &
!
                      ,advmnt_g(ngrid)%u3d  &
    		      ,advmnt_g(ngrid)%v3d  &
    		      ,advmnt_g(ngrid)%w3d  )
                                             
 
    if(theor_wind == on) then
 
        call prepare_theor_winds(dtlt,m1,m2,m3,ia,iz,ja,jz     &
                      ,advmnt_g(ngrid)%u3d  &
    		      ,advmnt_g(ngrid)%v3d  &
    		      ,advmnt_g(ngrid)%w3d  &
                      ,grid_g(ngrid)%dxt    &
                      ,grid_g(ngrid)%dyt    &
                      ,advmnt_g(ngrid)%dd0_3d  &
    		      ,advmnt_g(ngrid)%dd0_3du &
    		      ,advmnt_g(ngrid)%dd0_3dv &
		      ,advmnt_g(ngrid)%dd0_3dw )
    endif
    !- prepare Walcek's air densities
    
    call get_Walceks_densities(dtlt,m1,m2,m3 &
                      ,advmnt_g(ngrid)%u3d  &
    		      ,advmnt_g(ngrid)%v3d  &
    		      ,advmnt_g(ngrid)%w3d  &
                      ,advmnt_g(ngrid)%dd0_3d  &
    		      ,advmnt_g(ngrid)%dd0_3du &
    		      ,advmnt_g(ngrid)%dd0_3dv &
		      ,advmnt_g(ngrid)%dd0_3dw &
                      ,advmnt_g(ngrid)%den0_3d &
    		      ,advmnt_g(ngrid)%den1_3d &
    		      ,advmnt_g(ngrid)%den2_3d &
		      ,advmnt_g(ngrid)%den3_3d &
    		      ,advmnt_g(ngrid)%dxtW &
    		      ,advmnt_g(ngrid)%dytW &
		      ,advmnt_g(ngrid)%dztW &
                      ,grid_g(ngrid)%dxt    &
                      ,grid_g(ngrid)%dyt    )

         !call advect_mnt_stalone()



    CALL copyInitFields(ngrid,m1,m2,m3)
  	 
!
!    PRINT *,'InitialFieldsCopied: ',mynum;call flush(6)
    
!    WRITE(ctime,fmt='(I4.4)') int(time)

!LFR>      IF(bigdump==1) CALL dumpBorders('t'//ctime,'i',m1, m2, m3,mynum, &
!LFR>       nRec_i, procRecv_i, tagRecv_i, &
!LFR>       iaRecv_i, izRecv_i, jaRecv_i, jzRecv_i, &
!LFR>       bufRecvStart_i, bufRecvLength_i, bufRecvTotalLength_i, &
!LFR>       nSnd_i, procSend_i, tagSend_i, &
!LFR>       iaSend_i, izSend_i, jaSend_i, jzSend_i, &
!LFR>       bufSendStart_i, bufSendLength_i, bufSendTotalLength_i)	 

    CALL InitialFieldsUpdate(m1,m2,m3,newm2(ngrid),newm3(ngrid),ngrid,mynum, &
	   nRec_i, procRecv_i, tagRecv_i, &
	   iaRecv_i, izRecv_i, jaRecv_i, jzRecv_i, &
	   bufRecvStart_i, bufRecvLength_i, bufRecvTotalLength_i, &
	   nSnd_i, procSend_i, tagSend_i, &
	   iaSend_i, izSend_i, jaSend_i, jzSend_i, &
	   bufSendStart_i, bufSendLength_i, bufSendTotalLength_i)
    

!    PRINT *,'InitialFieldsUpdated i: ',mynum;call flush(6)

!LFR>      IF(bigdump==1) CALL dumpBorders('t'//ctime,'j',m1, m2, m3,mynum, &
!LFR>       nRec_j, procRecv_j, tagRecv_j, &
!LFR>       iaRecv_j, izRecv_j, jaRecv_j, jzRecv_j, &
!LFR>       bufRecvStart_j, bufRecvLength_j, bufRecvTotalLength_j, &
!LFR>       nSnd_j, procSend_j, tagSend_j, &
!LFR>       iaSend_j, izSend_j, jaSend_j, jzSend_j, &
!LFR>       bufSendStart_j, bufSendLength_j, bufSendTotalLength_j)	 

    CALL InitialFieldsUpdate(m1,m2,m3,newm2(ngrid),newm3(ngrid),ngrid,mynum, &
	   nRec_j, procRecv_j, tagRecv_j, &
	   iaRecv_j, izRecv_j, jaRecv_j, jzRecv_j, &
	   bufRecvStart_j, bufRecvLength_j, bufRecvTotalLength_j, &
	   nSnd_j, procSend_j, tagSend_j, &
	   iaSend_j, izSend_j, jaSend_j, jzSend_j, &
	   bufSendStart_j, bufSendLength_j, bufSendTotalLength_j)
	   
 !   PRINT *,'InitialFieldsUpdated j: ',mynum;call flush(6)


     !- ready to do advection, loop over all scalars
     do n=1,num_scalar(ngrid)     

!srf - somente para gases e aerossois
!     do n=num_scalar(ngrid) - NSPECIES_TRANSPORTED +1,num_scalar(ngrid)
!srf - somente para gases e aerossois
     !do n=num_scalar(ngrid) - NSPECIES_TRANSPORTED +1,num_scalar(ngrid)
        !print*,'nscalar1=',n
	!PRINT *,scalar_tab(n,ngrid)%name
	if (scalar_tab(n,ngrid)%name /= 'COP') cycle
        !print*,'nscalar2=',n
!srf - somente para gases e aerossois
!srf - somente para gases e aerossois
        !PRINT *,n; CALL flush(6)


 !srf-falta incluir controle para shaved_eta e aerossois com sedimentacao


      scalarp => scalar_tab(n,ngrid)%var_p
      scalart => scalar_tab(n,ngrid)%var_t

!for intel - descomente
      call atob(mxyzp,scalarp,advmnt_g(ngrid)%vc3d_jn) 
!LFR>       print *,'LFR  (.)(.) atob - passou ', mynum,n,ngrid,num_scalar(ngrid) ; call flush(6)

      call advect_mnt(ngrid,m1,m2,m3,ia,iz,ja,jz,dtlt &
!for pgi - descomente
!                     ,scalarp                 & ! initial mass mixing ratio (mmxr)
!for intel - descomente
 ,advmnt_g(ngrid)%vc3d_jn                     & ! initial mass mixing ratio (mmxr)
 ,advmnt_g(ngrid)%vc3d_out                    & ! final mmxr due advection
! ,advmntNew_g(ngrid)%vc3d_x			  &
! ,advmntNew_g(ngrid)%vc3d_y		     &
!-- 
		     ,advmntNew_g(ngrid)%u3d	 & ! 3D U wind
		    ,advmntNew_g(ngrid)%v3d	& ! 3D V wind
		    ,advmntNew_g(ngrid)%w3d	& ! 3D SIGMAZ wind
		     ,advmntNew_g(ngrid)%dd0_3d  & ! density 0
		     ,advmntNew_g(ngrid)%den0_3d & ! density D0; Eqs. 5  (JGR 2000)
		    ,advmntNew_g(ngrid)%den1_3d & ! density D1
		    ,advmntNew_g(ngrid)%den2_3d & ! density D2
		    ,advmntNew_g(ngrid)%den3_3d & ! density D3
		    ,advmntNew_g(ngrid)%dxtW	&
		    ,advmntNew_g(ngrid)%dytW	&
		    ,advmntNew_g(ngrid)%dztW	&

!		      ,advmnt_g(ngrid)%u3d     & ! 3D U wind
!		      ,advmnt_g(ngrid)%v3d     & ! 3D V wind
!		      ,advmnt_g(ngrid)%w3d     & ! 3D SIGMAZ wind
!		      ,advmnt_g(ngrid)%dd0_3d  & ! density 0
!		      ,advmnt_g(ngrid)%den0_3d & ! density D0; Eqs. 5  (JGR 2000)
!		      ,advmnt_g(ngrid)%den1_3d & ! density D1
!		      ,advmnt_g(ngrid)%den2_3d & ! density D2
!		      ,advmnt_g(ngrid)%den3_3d & ! density D3
!		      ,advmnt_g(ngrid)%dxtW    &
!		      ,advmnt_g(ngrid)%dytW    &
!		      ,advmnt_g(ngrid)%dztW    &

		     ,mynum,n                    &
		     ,nRec_i,nSnd_i,nRec_j,nSnd_j                &
                     ,bufSendTotalLength_i,bufSendTotalLength_j  &
	             ,bufRecvTotalLength_i,bufRecvTotalLength_j)
      		     
      call advtndc(m1,m2,m3,ia,iz,ja,jz        &
!                  ,scalarp  ,scratch%scr1 (1)  &
                  ,scalarp  ,advmnt_g(ngrid)%vc3d_out  &
                  ,scalart  ,dtlt,mynum        )

        iteste=iteste+1
     enddo
!     PRINT * ,'Deallocate Comm',mynum;call flush(6)  
!     CALL CommDeAlloc()
!stop 333
 
  END SUBROUTINE advmnt_driver
 !----------------------------------------------------
  
  SUBROUTINE initialize_advmnt(ngrids, mmzp,mmxp,mmyp)
   implicit none 
   INTEGER , INTENT(IN) :: ngrids
   INTEGER , INTENT(IN) :: mmxp(ngrids)
   INTEGER , INTENT(IN) :: mmyp(ngrids)
   INTEGER , INTENT(IN) :: mmzp(ngrids)
   
   INTEGER :: ng
   
   !maxz=maxval(mmzp(1:ngrids))
   !maxx=maxval(mmxp(1:ngrids))
   !maxy=maxval(mmyp(1:ngrids))
   
   IF(ALLOCATED(advmnt_g)) THEN 
      PRINT *,'Error in initialize_advmnt, sub: radvc_mnt: advmnt_g already allocated!'
      RETURN
   END IF
   ALLOCATE (advmnt_g(ngrids))
   ALLOCATE (advmntNew_g(ngrids))

   do ng=1,ngrids
      ALLOCATE(advmnt_g(ng)%u3d    (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%u3d=0.
      ALLOCATE(advmnt_g(ng)%v3d    (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%v3d=0.
      ALLOCATE(advmnt_g(ng)%w3d    (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%w3d=0.
      
      ALLOCATE(advmnt_g(ng)%dd0_3d (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3d =0.
      ALLOCATE(advmnt_g(ng)%dd0_3du(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3du=0.
      ALLOCATE(advmnt_g(ng)%dd0_3dv(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3dv=0.
      ALLOCATE(advmnt_g(ng)%dd0_3dw(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3dw=0.
      										      
      ALLOCATE(advmnt_g(ng)%den0_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den0_3d=0.
      ALLOCATE(advmnt_g(ng)%den1_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den1_3d=0.
      ALLOCATE(advmnt_g(ng)%den2_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den2_3d=0.
      ALLOCATE(advmnt_g(ng)%den3_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den3_3d=0.

      ALLOCATE(advmnt_g(ng)%dxtW(mmxp(ng),mmyp(ng))); advmnt_g(ng)%dxtW=0.
      ALLOCATE(advmnt_g(ng)%dytW(mmxp(ng),mmyp(ng))); advmnt_g(ng)%dytW=0.
      ALLOCATE(advmnt_g(ng)%dztW(mmzp(ng)))                  ; advmnt_g(ng)%dztW=0.


      ALLOCATE(advmnt_g(ng)%vc3d_jn  (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%vc3d_jn =0.
      ALLOCATE(advmnt_g(ng)%vc3d_out (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%vc3d_out=0.
      ALLOCATE(advmnt_g(ng)%vc3d_x   (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%vc3d_x  =0.
      ALLOCATE(advmnt_g(ng)%vc3d_y   (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%vc3d_y  =0.

!-------

      ALLOCATE(advmntNew_g(ng)%u3d    (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%u3d=0.
      ALLOCATE(advmntNew_g(ng)%v3d    (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%v3d=0.
      ALLOCATE(advmntNew_g(ng)%w3d    (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%w3d=0.
      
      ALLOCATE(advmntNew_g(ng)%dd0_3d (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%dd0_3d =0.
      ALLOCATE(advmntNew_g(ng)%dd0_3du(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%dd0_3du=0.
      ALLOCATE(advmntNew_g(ng)%dd0_3dv(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%dd0_3dv=0.
      ALLOCATE(advmntNew_g(ng)%dd0_3dw(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%dd0_3dw=0.
      										      
      ALLOCATE(advmntNew_g(ng)%den0_3d(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%den0_3d=0.
      ALLOCATE(advmntNew_g(ng)%den1_3d(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%den1_3d=0.
      ALLOCATE(advmntNew_g(ng)%den2_3d(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%den2_3d=0.
      ALLOCATE(advmntNew_g(ng)%den3_3d(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%den3_3d=0.

      ALLOCATE(advmntNew_g(ng)%dxtW(newM2(ng),newM3(ng))); advmntNew_g(ng)%dxtW=0.
      ALLOCATE(advmntNew_g(ng)%dytW(newM2(ng),newM3(ng))); advmntNew_g(ng)%dytW=0.
      ALLOCATE(advmntNew_g(ng)%dztW(mmzp(ng)))                  ; advmntNew_g(ng)%dztW=0.


      ALLOCATE(advmntNew_g(ng)%vc3d_jn  (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%vc3d_jn =0.
      ALLOCATE(advmntNew_g(ng)%vc3d_out (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%vc3d_out=0.
      ALLOCATE(advmntNew_g(ng)%vc3d_x   (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%vc3d_x  =0.
      ALLOCATE(advmntNew_g(ng)%vc3d_y   (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%vc3d_y  =0.

   enddo

   



  END SUBROUTINE initialize_advmnt
 !----------------------------------------------------
 
  SUBROUTINE initialize_densities(m1,m2,m3,dn0,dn0u,dn0v &
			    ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw )
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: m1,m2,m3
   REAL,DIMENSION(m1,m2,m3),intent(IN) ::dn0,dn0u,dn0v
   REAL,DIMENSION(m1,m2,m3),intent(OUT)::dd0_3d,dd0_3du,dd0_3dv,dd0_3dw
   ! local var
   integer i,j,k
   
   dd0_3d (:,:,:)=  dn0 (:,:,:) 
   dd0_3du(:,:,:)=  dn0u(:,:,:)
   dd0_3dv(:,:,:)=  dn0v(:,:,:)
   do j = 1,m3
      do i = 1,m2
         do k = 1,m1-1
            dd0_3dw(k,i,j) = 0.5*(dn0(k,i,j) +dn0(k+1,i,j))
         enddo
	 dd0_3dw(m1,i,j)=dd0_3dw(m1-1,i,j)
   enddo;enddo
   

  END SUBROUTINE initialize_densities
 !----------------------------------------------------
 
  SUBROUTINE initialize_grid_spacings(ng,m1,m2,m3,dxt,dyt,fmapt,rtgt &
			    ,dxtW,dytW,dztW )
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: ng,m1,m2,m3
   REAL,DIMENSION(m2,m3),intent(IN) :: dxt,dyt,fmapt,rtgt
   
   REAL,DIMENSION(m2,m3),intent(OUT):: dxtW,dytW
   REAL,DIMENSION(m1),intent(OUT):: dztW
   ! local var
   integer i,j,k
   real rtgti
   
   do j = 1,m3
      do i = 1,m2
   	rtgti = 1. / rtgt(i,j)
        
	!- at init/rams_grid.f90:
        !     dxt(i,j)=fmapt(i,j)/(xmn(i,ngrid)-xmn(i-1,ngrid))  
        !     dyt(i,j)=fmapt(i,j)/(ymn(j,ngrid)-ymn(j-1,ngrid))

   	 dxtW(i,j) = 1./(dxt(i,j) * fmapt(i,j) * rtgti)
   	 dytW(i,j) = 1./(dyt(i,j) * fmapt(i,j) * rtgti)
   enddo;enddo
   do k = 1,m1
    !- at init/gridset.f90:	
    !  dztn(k,ifm) = 1. / (zmn(k,ifm) - zmn(k-1,ifm))
    ! Por que o Jacobiano nao depende de Z, o dztw depende somente 
    ! de z.
    !dztW(k,i,j) = 1./ ( dzt(k) * rtgti * fmapt(i,j)**2 )
     dztW(k)	 = 1./ ( dztn(k,ng) ) !
 
   enddo
  END SUBROUTINE initialize_grid_spacings
 !----------------------------------------------------
 
  SUBROUTINE get_true_densities(m1,m2,m3,level,rtp,rv,pp,pi0,theta &
			       ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw )
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: m1,m2,m3,level
   REAL,DIMENSION(m1,m2,m3),intent(IN) ::rtp,rv,pp,pi0,theta
   REAL,DIMENSION(:,:,:),intent(OUT):: dd0_3d
   REAL,DIMENSION(m1,m2,m3),intent(OUT)  :: dd0_3du,dd0_3dv,dd0_3dw
   ! local var
   integer i,j,k
   real c3
   
   c3 = c2 * (cpi**c1)

   !- true air density at points "T"
   
   if( level == 0 ) then
     dd0_3d(:,:,:) = (c3/theta(:,:,:))*(pi0(:,:,:)+pp(:,:,:))**c1
   else
   do j = 1,m3
      do i = 1,m2
         do k = 1,m1
            dd0_3d(k,i,j) = (c3/theta(k,i,j))* (1. + rtp(k,i,j))/ &
	            (1. + 1.61*rv(k,i,j))*(pi0(k,i,j)+pp(k,i,j))**c1
         end do
      end do
   end do
   endif 

   !- true air density at points "U", "V" and "W":
   
   call fill_dn0uv(m1,m2,m3,dd0_3d,dd0_3du,dd0_3dv)

   do j = 1,m3
      do i = 1,m2
         do k = 1,m1-1
            dd0_3dw(k,i,j) = 0.5*(dd0_3d(k,i,j) + dd0_3d(k+1,i,j))
         enddo
	 dd0_3dw(m1,i,j)=dd0_3dw(m1-1,i,j)
   enddo;enddo

 END SUBROUTINE get_true_densities

 !----------------------------------------------------

 SUBROUTINE prepare_winds(dtlt,m1,m2,m3,ia,iz,ja,jz &
                            ,uc,up,vc,vp,wc,wp &
  		            ,fmapui &
		            ,fmapvi &
                            ,rtgt   &
                            ,rtgu   &
		            ,rtgv   &
		            ,f13t   &
                            ,f23t   &
			    ,u3d,v3d,w3d)
   
   implicit none 
   INTEGER , INTENT(IN) :: m1,m2,m3,ia,iz,ja,jz
   REAL, INTENT(IN) :: dtlt
   REAL,DIMENSION(m1,m2,m3),intent(INOUT) :: uc,up,vc,vp,wc,wp
!   REAL,DIMENSION(m1,m2,m3),intent(IN) :: uc,up,vc,vp,wc,wp ! TMP <<<<<<<<<<<<<<<<<<<<<<<
   REAL,DIMENSION(m2,m3)   ,intent(IN) :: rtgt,rtgu,rtgv,fmapui,fmapvi,f13t,f23t
   
   REAL,DIMENSION(m1,m2,m3),intent(OUT)::u3d,v3d,w3d
  
   !- local var
   real   dtlto2
   integer jm,jp,im,ip 
   integer i,j,k
   real :: c1,c2,rtgti

   dtlto2 = .5! * dtlt

  ! u3d, u3d, and w3d are input as the velocity components (averaged
  ! between past and current time levels) times dtlt.
   do j=1,m3
     do i = 1,m2
      do k = 1,m1
  
          w3d(k,i,j) = ( wc(k,i,j) + wp(k,i,j) )*dtlto2
          u3d(k,i,j) = ( uc(k,i,j) + up(k,i,j) )*dtlto2
          v3d(k,i,j) = ( vc(k,i,j) + vp(k,i,j) )*dtlto2
  
   enddo;enddo;enddo


  !return ! for pure cartesian coordinates

  ! here w3d is the cartesian vertical velocity 
  
  ! Add contribution to w3d from horiz winds crossing sloping sigma surfaces,
  ! and include 1/rtgt factor in w3d
  do j = 1,m3		 
     jm = max(1,j-1)
     jp = min(m3,j+1)
     do i = 1,m2
  	im = max(1,i-1)
  	ip = min(m2,i+1)
        rtgti = 1. / rtgt(i,j)
	
	do k = 1,m1-1
	    w3d(k,i,j) = ( (u3d(k,i,j) + u3d(k+1,i,j) + u3d(k,im,j) + u3d(k+1,im,j) ) * f13t(i,j)  &
	               +   (v3d(k,i,j) + v3d(k+1,i,j) + v3d(k,i,jm) + v3d(k+1,i,jm) ) * f23t(i,j)  &
		         ) * hw4(k)  &
	               + w3d(k,i,j) * rtgti

  	enddo
     enddo
  enddo
  ! here w3d is the sigma_z velocity 

  !- including map factors on U,V:
  do j = 1,m3
     do i = 1,m2
  	c1 = fmapui(i,j) * rtgu(i,j)
  	c2 = fmapvi(i,j) * rtgv(i,j)
  	do k = 1,m1-1
  	   u3d(k,i,j) = u3d(k,i,j) * c1 
  	   v3d(k,i,j) = v3d(k,i,j) * c2 
	   
	   
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	   !u3d(k,i,j) = 0. !<<<<<<<<<<<<<<<<<<<<<<
	   
	   
	   !v3d(k,i,j) = -1. !<<<<<<<<<<<<<<<<<<<<<<

           !w3d(k,i,j) =	 0.
	   !wc(k,i,j)=w3d(k,i,j)
	   !uc(k,i,j)=u3d(k,i,j)
	   !vc(k,i,j)=v3d(k,i,j)
	   !wp(k,i,j)=w3d(k,i,j)
	   !up(k,i,j)=u3d(k,i,j)
	   !vp(k,i,j)=v3d(k,i,j)
	   
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  	enddo
     enddo
  enddo

  END SUBROUTINE prepare_winds
 !----------------------------------------------------

  SUBROUTINE get_Walceks_densities(dt,m1,m2,m3,u3d,v3d,w3d &
			    ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw &
			    ,den0_3d,den1_3d,den2_3d,den3_3d &
			    ,dxtW,dytW,dztW,dxt,dyt)

   IMPLICIT NONE 
   !-in
   INTEGER , INTENT(IN) :: m1,m2,m3
   REAL    , INTENT(IN) :: dt
   REAL, DIMENSION(m1)      , INTENT(IN) :: dztW
   REAL, DIMENSION(m2,m3)   , INTENT(IN) :: dxtW,dytW

   REAL, DIMENSION(m2,m3)   , INTENT(IN) :: dxt,dyt

   REAL, DIMENSION(m1,m2,m3), INTENT(IN) :: u3d,v3d,w3d     &
			                   ,dd0_3du &
					   ,dd0_3dv,dd0_3dw
   REAL, DIMENSION(m1,m2,m3), INTENT(INOUT) :: dd0_3d 
   !-out
   REAL, DIMENSION(m1,m2,m3), INTENT(OUT):: den0_3d,den1_3d &
                                           ,den2_3d,den3_3d
   
   ! local var
   integer i,j,k
   
    DO  j=m3,2,-1
      DO  i=2,m2
        DO k = 2,m1
	

!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!dd0_3d(k,i,j)=1.
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<


            den0_3d(k,i,j)=dd0_3d(k,i,j)



	    
!	    den1_3d(k,i,j)=den0_3d(k,i,j) - dt/dxt(i,j)*(d0         *u3d(k,i,j)- d0           *u3d(k,i-1,j ))
!	    den1_3d(k,i,j)=den0_3d(k,i,j) - dt*dxt(i,j)*(dd0_3du(k,i,j)*u3d(k,i,j)- dd0_3du(k,i-1,j)*u3d(k,i-1,j  ))
	    den1_3d(k,i,j)=den0_3d(k,i,j) - dt/dxtW(i,j)*(dd0_3du(k,i,j)*u3d(k,i,j)- dd0_3du(k,i-1,j)*u3d(k,i-1,j  ))

!           den2_3d(k,i,j)=den1_3d(k,i,j)- dt/dx*(d0*v3d(k,i,j)- d0*v3d(k,i,j-1 ))
!           den2_3d(k,i,j)=den1_3d(k,i,j)- dt*dyt(i,j)*(dd0_3dv(k,i,j)*v3d(k,i,j)- dd0_3dv(k,i,j-1)*v3d(k,i,j-1  ))
            den2_3d(k,i,j)=den1_3d(k,i,j)- dt/dytW(i,j)*(dd0_3dv(k,i,j)*v3d(k,i,j)- dd0_3dv(k,i,j-1)*v3d(k,i,j-1  ))
!
!           den3_3d(k,i,j)=den2_3d(k,i,j)-dt/dx         *(d0*w3d(k,i,j)-d0*w3d(  k-1,i,j))
!           den3_3d(k,i,j)=den2_3d(k,i,j)-dt*dzt(k)     *(dd0_3dw(k,i,j)*w3d(k,i,j)-dd0_3dw(k-1,i,j)*w3d(k-1,i,j))
            den3_3d(k,i,j)=den2_3d(k,i,j)-dt/dztW(k)*(dd0_3dw(k,i,j)*w3d(k,i,j)-dd0_3dw(k-1,i,j)*w3d(k-1,i,j))

       END DO !i
      END DO !j
    END DO !k  

 

 END SUBROUTINE get_Walceks_densities
 !----------------------------------------------------
 SUBROUTINE advect_mnt(ngrid,m1,m2,m3,ia,iz,ja,jz,dt,vc3d_in,vc3d_out &  
                      ,u3d,v3d,w3d,dd0,den0,den1,den2,den3 &
		      ,dxtW,dytW,dztW,mynum,n, &
		      nRec_i,nSnd_i,nRec_j,nSnd_j, &
                      bufSendTotalLength_i,bufSendTotalLength_j, &
			bufRecvTotalLength_i,bufRecvTotalLength_j)
   
   use node_mod, only : nodei0,nodej0
   use mem_grid, only : time ! tmp
  IMPLICIT NONE
  INTEGER , INTENT(IN) :: m1,ngrid
  INTEGER , INTENT(IN) :: m2
  INTEGER , INTENT(IN) :: m3
  INTEGER , INTENT(IN) :: ia
  INTEGER , INTENT(IN) :: iz
  INTEGER , INTENT(IN) :: ja
  INTEGER , INTENT(IN) :: jz,n
  INTEGER , INTENT(IN) :: mynum
  INTEGER, INTENT(IN)  :: nRec_i,nSnd_i,nRec_j,nSnd_j
  INTEGER, INTENT(IN)  :: bufSendTotalLength_i,bufSendTotalLength_j
  INTEGER, INTENT(IN)  :: bufRecvTotalLength_i,bufRecvTotalLength_j
  REAL    , INTENT(IN) :: dt
  
  !REAL,DIMENSION(m1)      , INTENT(IN)  :: dzt
  !REAL,DIMENSION(m2,m3)   , INTENT(IN)  :: dxt,dyt,rtgt
  
  REAL,DIMENSION(m1),       INTENT(IN) :: dztW 
  REAL,DIMENSION(newM2(ngrid),newM3(ngrid))   , INTENT(IN) :: dxtW,dytW

  REAL,DIMENSION(m1,m2,m3), INTENT(IN) :: vc3d_in! initial mixing ratio 
  REAL,DIMENSION(m1,newM2(ngrid),newM3(ngrid)), INTENT(IN) :: u3d,v3d,w3d
  REAL,DIMENSION(m1,newM2(ngrid),newM3(ngrid)), intent(IN) :: dd0,den0,den1,den2,den3

  !REAL,DIMENSION(m1,m2,m3), INTENT(OUT) :: vc3d_x,vc3d_y
  REAL,DIMENSION(m1,m2,m3), INTENT(OUT) :: vc3d_out ! updated mixing ratio due advection

  !- local var
  REAL,DIMENSION(m1)               :: dxx  
  REAL,DIMENSION(m2,m3)            :: dxy  
  real masscon,initialmass,vol
  integer :: i,j,k,nrec,it
  
  CHARACTER(LEN=2) :: nf

  REAL,DIMENSION(m1,m2,m3) :: vc3dout_X,vc3dout_Y,vc3dout_Z
  REAL,DIMENSION(m1,newM2(ngrid),newM3(ngrid)) :: vc3X,vc3Y,vc3Z  

  INTEGER :: iBegin, iEnd, jBegin, jEnd
     
     iBegin=newIa(ngrid)-1
     iEnd=iBegin+m2-1
     jBegin=newJa(ngrid)-1
     jEnd=jBegin+m3-1
          

      !WRITE (71,*) m2,m3,newM2(ngrid),newM3(ngrid),newIa(ngrid),newJa(ngrid)
      !WRITE (71,*) iBegin,iEnd,jBegin,jEnd
      !CALL flush(71)

      vc3z(1:m1,iBegin:iEnd,jBegin:jEnd)=vc3d_in(1:m1,1:m2,1:m3)
!    vc3z=float(mynum)
    
! write(60+mynum,*) '1', vc3z(4,1:newM2(ngrid),7); call flush(60+mynum)
    
     call UpdateBorders(m1, newm2(ngrid), newm3(ngrid), vc3z, &
          nRec_j, procRecv_j, tagRecv_j, &
          iaRecv_j, izRecv_j, jaRecv_j, jzRecv_j, &
          bufRecvStart_j, bufRecvLength_j, bufRecvTotalLength_j, &
          nSnd_j, procSend_j, tagSend_j, &
          iaSend_j, izSend_j, jaSend_j, jzSend_j, &
          bufSendStart_j, bufSendLength_j, bufSendTotalLength_j)
!write(60+mynum,*) '2',vc3z(4,1:newM2(ngrid),7); call flush(60+mynum)

     call Advec3d_X(m1,newM2(ngrid),newM3(ngrid),2,newM2(ngrid)-1,2,newM3(ngrid)-1,vc3z,u3d,den0, &
                 den1,dt,dxtW,dd0,vc3x,mynum )
		 
!     call Advec3d_X(m1, newm2(ngrid), newm3(ngrid), &
!          vc3z, u3d, den0, den1, dt, dxtW, dd0, vc3x)
!srf tmp
!     vc3d_out(1:m1,1:m2,1:m3)=vc3x(1:m1,iBegin:iEnd,jBegin:jEnd)
!write(60+mynum,*) vc3x(4,1:newM2(ngrid),10); call flush(60+mynum)
 !    return
!srf tmp


     ! Do J-advection next over all I-k lines
     call UpdateBorders(m1, newm2(ngrid), newm3(ngrid), vc3x, &
	 nRec_i, procRecv_i, tagRecv_i, &
	 iaRecv_i, izRecv_i, jaRecv_i, jzRecv_i, &
	 bufRecvStart_i, bufRecvLength_i, bufRecvTotalLength_i, &
	 nSnd_i, procSend_i, tagSend_i, &
	 iaSend_i, izSend_i, jaSend_i, jzSend_i, &
	 bufSendStart_i, bufSendLength_i, bufSendTotalLength_i)

    call Advec3d_Y(m1,newM2(ngrid),newM3(ngrid),2,newM2(ngrid)-1,2,newM3(ngrid)-1,vc3x ,v3d,den1, &
                 den2,dt,dytW,dd0,vc3y )


!srf tmp
     vc3d_out(1:m1,1:m2,1:m3)=vc3y(1:m1,iBegin:iEnd,jBegin:jEnd)
     return
!srf tmp


!     call Advec3d_Y(m1, newm2(ngrid), newm3(ngrid), &
!          vc3x, v3d, den1, den2, dt, dytW, dd0, vc3y)

        !--- do k-advection next over all I-J lines
        ! vc3dout_Z = razao de mistura final (apos a adveccao nas 3 dimensoes)
             call Advec3d_Z(m1,newM2(ngrid),newM3(ngrid),2,newM2(ngrid)-1,2,newM3(ngrid)-1,vc3y, &
	                 w3d,den2,den3,dt,dztW,dd0,vc3z)
!LFR>         call Advec3d_Z(m1, newm2(ngrid), newm3(ngrid), &
!LFR>              vc3y, w3d, den2, den3, dt, dztW, dd0, vc3z)

   
     vc3d_out(1:m1,1:m2,1:m3)=vc3z(1:m1,iBegin:iEnd,jBegin:jEnd)

 END SUBROUTINE advect_mnt
 !----------------------------------------------------
 
 
  SUBROUTINE prepare_theor_winds(dtlt,m1,m2,m3,ia,iz,ja,jz &
                            ,u3d,v3d,w3d&
			    ,dxt,dyt &
			    ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw )

   
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: m1,m2,m3,ia,iz,ja,jz
   REAL, INTENT(IN) :: dtlt
   REAL, DIMENSION(m2,m3)   , INTENT(IN) :: dxt,dyt

   REAL,DIMENSION(m1,m2,m3),intent(OUT)::u3d,v3d,w3d     &
                                        ,dd0_3d ,dd0_3du &
					,dd0_3dv,dd0_3dw 
  
   !- local var
   REAL   :: dtlto2
   INTEGER :: i,j,k
   REAL  :: ai0s  =  25.0
   REAL  :: aj0s  =  50.0
   REAL  :: umx   =   10.0
   REAL,PARAMETER :: pii   =   3.141592653589793
   REAL    :: umax  =   0.0
   REAL    :: anrev,curnt,rx,xa,ilop,iwndty,nrec
   
   dtlto2 =  10.!*dtlt
    
     !WRITE(6,*) ' Wind fields?  0-rotating; or  1-divergent winds'
     iwndty = 0  ! 0-rotating
     !iwndty = 1  ! 1-divergent winds

     IF(iwndty==1) ai0s= 50.5
     ilop= ai0s-21.  ! needed for printouts
     ! Define wind fields (rotation or divergent) and initial mixing ratios
     !  Cone at (25,50) for rotating winds; Cone at (50,50) divergent winds
     DO k = 1,m1
   	DO  j=m3,1,-1
   	   DO  i=1,m2
   	   
   	       dd0_3d (k,i,j)=1.
	       dd0_3du(k,i,j)=1. 
               dd0_3dv(k,i,j)=1.
	       dd0_3dw(k,i,j)=1. 

   	       u3d(k,i,j)= -2.*umx*(REAL(j)-REAL(110)/2.-.5)/REAL(110)*dtlto2
   	       v3d(k,i,j)=  2.*umx*(REAL(i)-REAL(100)/2.-.5)/REAL(100)*dtlto2
 	       w3d(k,i,j)= 0.                                        *dtlto2
   	   
   	      IF(iwndty==1) THEN
   		 xa=pii/25.
   		 IF(J>0) u3d(k,i,j)=umx*SIN(xa*REAL(i))*SIN(xa*(REAL(j)))
   		 IF(I>0) v3d(k,i,j)=umx*COS(xa*(REAL(i)-.5))*COS(XA*(REAL(j)+.5))
   	      END IF
   	   
   	   
   	      umax= MAX(ABS(u3d(k,i,j)),ABS(v3d(k,i,j)),umax)
   	      rx= SQRT((REAL(i)-ai0s)**2.+(REAL(j)-aj0s)**2.)

   	   END DO !i
   	END DO !j
     END DO !k  

  END SUBROUTINE prepare_theor_winds
  
   subroutine UpdateBorders(m1, m2, m3, field, &
	nRecv, procRecv, tagRecv, iaRecv, izRecv, jaRecv, jzRecv, &
	bufRecvStart, bufRecvLength, bufRecvTotalLength, &
	nSend, procSend, tagSend, iaSend, izSend, jaSend, jzSend, &
	bufSendStart, bufSendLength, bufSendTotalLength)
     integer, intent(in) :: m1
     integer, intent(in) :: m2
     integer, intent(in) :: m3
     real,    intent(inout) :: field(m1,m2,m3)
     integer, intent(in) :: nRecv
     integer, intent(in) :: procRecv(nRecv)
     integer, intent(in) :: tagRecv(nRecv)
     integer, intent(in) :: iaRecv(nRecv)
     integer, intent(in) :: izRecv(nRecv)
     integer, intent(in) :: jaRecv(nRecv)
     integer, intent(in) :: jzRecv(nRecv)
     integer, intent(in) :: bufRecvStart(nRecv)
     integer, intent(in) :: bufRecvLength(nRecv)
     integer, intent(in) :: bufRecvTotalLength
     integer, intent(in) :: nSend
     integer, intent(in) :: procSend(nSend)
     integer, intent(in) :: tagSend(nSend)
     integer, intent(in) :: iaSend(nSend)
     integer, intent(in) :: izSend(nSend)
     integer, intent(in) :: jaSend(nSend)
     integer, intent(in) :: jzSend(nSend)
     integer, intent(in) :: bufSendStart(nSend)
     integer, intent(in) :: bufSendLength(nSend)
     integer, intent(in) :: bufSendTotalLength
 
 
     integer :: i, j, iCnt, i1, i2, cnt, iRecv, iSend, ierr
     real, allocatable :: bufRecv(:)
     integer :: reqRecv(nRecv)
     real, allocatable :: bufSend(:)
     integer :: reqSend(nSend)
     integer :: status(MPI_STATUS_SIZE)
     integer :: status2(MPI_STATUS_SIZE,nSend)
 
     allocate(bufRecv(bufRecvTotalLength))
     do iRecv = 1, nRecv
	call MPI_Irecv(bufRecv(bufRecvStart(iRecv)), bufRecvLength(iRecv), &
	     MPI_REAL, &
	     procRecv(iRecv), tagRecv(iRecv), MPI_COMM_WORLD, &
	     reqRecv(iRecv), ierr)
     end do
 
     allocate(bufSend(bufSendTotalLength))
     do iSend = 1, nSend
	i1 = bufSendStart(iSend)
	iCnt = bufSendStart(iSend)
	i2 = bufSendLength(iSend)
	do j = jaSend(iSend), jzSend(iSend)
	   do i = iaSend(iSend), izSend(iSend)
	      bufSend(iCnt:iCnt+m1-1) = field(1:m1,i,j)
	      iCnt = iCnt+m1
	   end do
	end do
	call MPI_Isend(bufSend(i1), i2, MPI_REAL, &
	     procSend(iSend), tagSend(iSend), MPI_COMM_WORLD, &
	     reqSend(iSend), ierr)
     end do
 
     do cnt = 1, nRecv
	call MPI_Waitany(nRecv, reqRecv, iRecv, status, ierr)
	i1 = bufRecvStart(iRecv)
	iCnt = bufRecvStart(iRecv)
	i2 = bufRecvLength(iRecv)
	do j = jaRecv(iRecv), jzRecv(iRecv)
	   do i = iaRecv(iRecv), izRecv(iRecv)
	      field(1:m1,i,j) = bufRecv(iCnt:iCnt+m1-1)
	      iCnt = iCnt+m1
	   end do
	end do
     end do
 
     call MPI_Waitall(nSend, reqSend, status2, ierr)
 
     deallocate(bufSend)
     deallocate(bufRecv)
   end subroutine UpdateBorders

   SUBROUTINE InitCommMatrix(Mynum,master_num)
   
      use rpara, only: &
               nmachs,ixb,ixe,iyb,iye
      use mem_grid, only: &
               ngrids, nnxp,nnyp
	       
      IMPLICIT NONE
      
      INCLUDE 'mpif.h'
  
      INTEGER, INTENT(IN) :: mynum,master_num
      
      INTEGER :: i,j,ngrid,nproc,nprocrec,fromto,ncount,nvar,ierr,inx,iny,ng
      INTEGER :: procmap1d(ngrids*maxval(nnxp)*maxval(nnyp))
      INTEGER :: info(1),ihandle(2),countItag
      
      INTEGER  :: tagcount,jj
      INTEGER,DIMENSION(ngrids) :: nRecv_i,nRecv_j
      INTEGER,DIMENSION(ngrids,maxval(nnxp)*maxval(nnyp)) :: procRecv_i,procSend_i
      INTEGER,DIMENSION(ngrids,maxval(nnxp)*maxval(nnyp)) :: procRecv_j,procSend_j
      INTEGER,DIMENSION(ngrids,maxval(nnxp)*maxval(nnyp)) :: tagSend_i,tagSend_j
      INTEGER,DIMENSION(ngrids,maxval(nnxp)*maxval(nnyp)) :: tagRecv_i,tagRecv_j
 
      INTEGER :: unitout
      CHARACTER(LEN=2) :: c_mynum
 
      ALLOCATE(procMap(ngrids,maxval(nnxp),maxval(nnyp)));procMap=0
      ALLOCATE(handle(nmachs+100000))      
      ALLOCATE(west(ngrids),east(ngrids),north(ngrids),south(ngrids))
      IF(mynum==master_num) THEN !I am the master
        
         PRINT *,'Notice from the master:'
         PRINT *,'Number of processors: ',nmachs
         PRINT *,'Master node is the  : ',mynum
         PRINT *,'ngrids, maxval(nnxp),maxval(nnyp)=',ngrids, maxval(nnxp),maxval(nnyp)

     
	 info(1)=nmachs
         CALL flush(6)    

         DO ngrid=1,ngrids
            DO nproc=1,nmachs
               DO i=ixb(nproc,ngrid),ixe(nproc,ngrid)
                  DO j=iyb(nproc,ngrid),iye(nproc,ngrid)
	             procMap(ngrid,i,j)=nproc   
	          END DO
               END DO
	    END DO
         END DO
	
	  WRITE(*,FMT='(A)') '--(.)(.)-- Map of processors for advection --(.)(.)--';CALL flush(6)
	  DO i=1,nnxp(1)
	     WRITE(*,FMT='(A,I3.3,2X,$)') 'I= ',i
	     DO j=1,nnyp(1)
		WRITE(*,FMT='(I2.2,A,$)') procMap(1,i,j),' ';CALL flush(6)
	     END DO
	     WRITE(*,FMT='(A)') '';CALL flush(6)
	  END DO


         ncount=0
	 DO ngrid=1,ngrids
	    DO i=1,maxval(nnxp)
	       DO j=1,maxval(nnyp)
	             ncount=ncount+1
	             procmap1d(ncount)=procMap(ngrid,i,j)
	       END DO
	    END DO
	 END DO
	 DO nproc=1,nmachs	
	    !PRINT *,'Notice from ',mynum,' procmap sended to proc #',nproc; CALL flush(6)    
	    CALL send(info,1,nproc,nproc+50000)
	    CALL send(procmap1d,ngrids*maxval(nnxp)*maxval(nnyp),nproc,nproc)
	 END DO


      ELSE !I am a slave      


	 CALL get(info,1,0,mynum+50000)
	 CALL get(procmap1d,ngrids*maxval(nnxp)*maxval(nnyp),master_num,mynum)
	 CALL parf_wait_nostatus(handle(mynum+50000))
	 CALL parf_wait_nostatus(handle(mynum))

         ncount=0
         DO ngrid=1,ngrids
	    DO i=1,maxval(nnxp)
	       DO j=1,maxval(nnyp)
	             ncount=ncount+1
	             procMap(ngrid,i,j)=procmap1d(ncount)
	          END DO
	    END DO
         END DO      
      END IF

      CALL PrepareComm(procMap,mynum)
      
      IF(bigdump==1) THEN 
         unitout=20+mynum
	 write(c_mynum,FMT='(I2.2)') mynum 
         OPEN (UNIT=unitout,file='commDump'//c_mynum//'.dat')
         DO ng=1,ngrids
            WRITE (unitout,FMT='(A)') '==== i direction ====' 
	    WRITE (unitout,FMT='(A5,1X,11(A2,1X))') &
	    ' Tag ','sp','rp','ln','ia','iz','ja','jz','ia','iz','ja','jz'
	    CALL flush(unitout)
            DO i=1,nSend_i
	       WRITE (unitout,FMT='(I5.5,1X,11(I2.2,1X))') &
	       comm_i(ng)%tag(i), &
	       comm_i(ng)%procSend(i), &
	       comm_i(ng)%procRecv(i), &
	       comm_i(ng)%length(i), &
	       comm_i(ng)%iaSend(i), &
	       comm_i(ng)%izSend(i), &
	       comm_i(ng)%jaSend(i), &
	       comm_i(ng)%jzSend(i), &
	       comm_i(ng)%iaRecv(i), &
	       comm_i(ng)%izRecv(i), &
	       comm_i(ng)%jaRecv(i), &	 
	       comm_i(ng)%jzRecv(i)
	       CALL flush(unitout)    
	    END DO
         END DO   
         DO ng=1,ngrids
            WRITE (unitout,FMT='(A)') '==== j direction ====' 
	    WRITE (unitout,FMT='(A5,1X,11(A2,1X))') &
	    ' Tag ','sp','rp','ln','ia','iz','ja','jz','ia','iz','ja','jz'
	    CALL flush(unitout)
            DO i=1,nSend_j
	       WRITE (unitout,FMT='(I5.5,1X,11(I2.2,1X))') &
	       comm_j(ng)%tag(i), &
	       comm_j(ng)%procSend(i), &
	       comm_j(ng)%procRecv(i), &
	       comm_j(ng)%length(i), &
	       comm_j(ng)%iaSend(i), &
	       comm_j(ng)%izSend(i), &
	       comm_j(ng)%jaSend(i), &
	       comm_j(ng)%jzSend(i), &
	       comm_j(ng)%iaRecv(i), &
	       comm_j(ng)%izRecv(i), &
	       comm_j(ng)%jaRecv(i), &   
	       comm_j(ng)%jzRecv(i) 
	       CALL flush(unitout)   
	    END DO	    
         END DO 
         DO ng=1,ngrids	 
	    WRITE (unitout,FMT='("Neibour (Grid,W,E,N,S):",4(L1,1X))') ng,west(ng),east(ng),north(ng),south(ng)
         END DO
	 CLOSE(UNIT=unitout)  
      END IF
      
     
   END SUBROUTINE InitCommMatrix  
  !-------------------------------------------------------

  SUBROUTINE get_real(var,m1,from,tag)
       
     INTEGER, INTENT(IN) :: from,tag
     INTEGER,INTENT(IN)  :: m1
     REAL,INTENT(OUT)    :: var(m1)
     INTEGER :: Ireq
     INCLUDE "i8.h"
     
     CALL parf_get_noblock_real(var, &
	                   m1, &
	                   from, &
	                   tag, &
			   handle(tag)) 
  !call start_watch(w(1))     ! This starts the watch
  !   CALL parf_wait_nostatus(handle(tag))
  !call stop_watch(w(1))      ! This stops the watch
     
  END SUBROUTINE get_real
  
  SUBROUTINE get_int(var,m1,from,tag)
       
     INTEGER, INTENT(IN) :: from,tag
     INTEGER,INTENT(IN)  :: m1
     INTEGER,INTENT(OUT)    :: var(m1)
     INTEGER :: Ireq
     INCLUDE "i8.h"
     
     CALL parf_get_noblock_int(var, &
	                   m1, &
	                   from, &
	                   tag, &
			   handle(tag)) 
  !call start_watch(w(1))     ! This starts the watch
  !   CALL parf_wait_nostatus(handle(tag))
  !call stop_watch(w(1))      ! This stops the watch
     
  END SUBROUTINE get_int

  SUBROUTINE get_logical(var,m1,from,tag)

     INTEGER, INTENT(IN) :: from,tag
     INTEGER,INTENT(IN)  :: m1
     LOGICAL,INTENT(OUT) :: var(m1)
     INTEGER :: Ireq
     INCLUDE "i8.h"
     
     CALL parf_get_noblock_logical(var, &
	                   m1, &
	                   from, &
	                   tag, &
			   handle(tag)) 
  !call start_watch(w(1))     ! This starts the watch
  !   CALL parf_wait_nostatus(handle(tag))
  !call stop_watch(w(1))      ! This stops the watch
     
  END SUBROUTINE get_logical
  
  SUBROUTINE send_real(var,m1,from,tag)

     INTEGER, INTENT(IN) :: from,tag
     INTEGER,INTENT(IN)  :: m1
     REAL,INTENT(IN)    :: var(m1)
     INTEGER :: Ireq
     INCLUDE "i8.h"

     CALL parf_send_noblock_real(var, &
	                  m1, &
	                  from, &
			  tag, &
			  handle(tag)) 

  END SUBROUTINE send_real  
  
  SUBROUTINE send_int(var,m1,from,tag)

     INTEGER, INTENT(IN) :: from,tag
     INTEGER,INTENT(IN)  :: m1
     INTEGER,INTENT(IN)    :: var(m1)
     INTEGER :: Ireq
     INCLUDE "i8.h"

     CALL parf_send_noblock_int(var, &
	                  m1, &
	                  from, &
			  tag, &
			  handle(tag)) 

  END SUBROUTINE send_int 

  SUBROUTINE send_logical(var,m1,from,tag)

     INTEGER, INTENT(IN) :: from,tag
     INTEGER,INTENT(IN)  :: m1
     LOGICAL,INTENT(IN)    :: var(m1)
     INTEGER :: Ireq
     INCLUDE "i8.h"

     CALL parf_send_noblock_logical(var, &
	                  m1, &
	                  from, &
			  tag, &
			  handle(tag)) 

  END SUBROUTINE send_logical 
  

  SUBROUTINE PrepareComm(procMap,mynum)
      use mem_grid, only: &
               ngrids, nnxp,nnyp

      use node_mod, only : nodei0,nodej0

    INTEGER,INTENT(IN) :: mynum
    INTEGER,INTENT(IN),DIMENSION(ngrids,maxval(nnxp),maxval(nnyp)) :: procmap
    
    INTEGER :: ng,i,ii,j,jj,p1,p2,ini,fin,tagCount
    INTEGER :: aux,nmess_i,nmess_j,istep,jstep,igz
    
    
    nMess_i=0
    nMess_j=0
    tagCount=0
    west=.false.
    east=.false.
    south=.false.
    north=.false.
    
    ALLOCATE(comm_i(ngrids),comm_j(ngrids))
    
    !Counting the numbers of communications i-i and j-j
    DO ng=1,ngrids
       i=1
       DO WHILE (i<=nnxp(ng))
          i=i+1
          DO j=2,nnyp(ng)
             !Se existe diferenca entre esse ponto e o anterior
	     IF (procMap(ng,i,j)/=procMap(ng,i,j-1) .and. procMap(ng,i,j-1)/=0 &
	         .and. procMap(ng,i,j)/=0) THEN
                p1=procMap(ng,i,j)
		p2=procMap(ng,i,j-1)
		ini=i
		!existe comunicacao entre:
                !p1 -> p2
                !p2 -> p1		
		IF(procMap(ng,i,j)==myNum .or. procMap(ng,i,j-1)==myNum) &
		       nMess_i=nMess_i+1
		!verificando ateh onde p1 fala com p2
		DO ii=i+1,nnxp(ng)
		   !Se um dos dois eh diferente acabou
		   IF(procMap(ng,ii,j)/=p1 .or. procMap(ng,ii,j-1)/=p2) THEN
		      i=ii-1
		      EXIT
		   END IF
		   fin=ii
		END DO 
             END IF
	  END DO
       END DO
       ALLOCATE (comm_i(ng)%tag(nMess_i*2))
       ALLOCATE (comm_i(ng)%length(nMess_i*2))
       ALLOCATE (comm_i(ng)%procSend(nMess_i*2))
       ALLOCATE (comm_i(ng)%iaSend(nMess_i*2))
       ALLOCATE (comm_i(ng)%izSend(nMess_i*2))
       ALLOCATE (comm_i(ng)%jaSend(nMess_i*2))
       ALLOCATE (comm_i(ng)%jzSend(nMess_i*2))
       ALLOCATE (comm_i(ng)%procRecv(nMess_i*2))
       ALLOCATE (comm_i(ng)%iaRecv(nMess_i*2))
       ALLOCATE (comm_i(ng)%izRecv(nMess_i*2))
       ALLOCATE (comm_i(ng)%jaRecv(nMess_i*2))
       ALLOCATE (comm_i(ng)%jzRecv(nMess_i*2))
    END DO
    DO ng=1,ngrids
       j=1
       DO WHILE (j<=nnyp(ng))
          j=j+1
          DO i=2,nnxp(ng)
             !Se existe diferenca entre esse ponto e o anterior
	     IF (procMap(ng,i,j)/=procMap(ng,i-1,j) .and. procMap(ng,i-1,j)/=0 &
	         .and. procMap(ng,i,j)/=0) THEN
                p1=procMap(ng,i,j)
		p2=procMap(ng,i-1,j)
		ini=j
		!existe comunicacao entre:
                !p1 -> p2
                !p2 -> p1		
		IF(procMap(ng,i,j)==myNum .or. procMap(ng,i-1,j)==myNum) &
		       nMess_j=nMess_j+1
		!verificando ateh onde p1 fala com p2
		DO jj=j+1,nnyp(ng)
		   !Se um dos dois eh diferente acabou
		   IF(procMap(ng,i,jj)/=p1 .or. procMap(ng,i-1,jj)/=p2) THEN
		      j=jj-1
		      EXIT
		   END IF
		   fin=jj
		END DO 
             END IF
	  END DO
       END DO
       ALLOCATE (comm_j(ng)%tag(nMess_j*2))
       ALLOCATE (comm_j(ng)%length(nMess_j*2))
       ALLOCATE (comm_j(ng)%procSend(nMess_j*2))
       ALLOCATE (comm_j(ng)%iaSend(nMess_j*2))
       ALLOCATE (comm_j(ng)%izSend(nMess_j*2))
       ALLOCATE (comm_j(ng)%jaSend(nMess_j*2))
       ALLOCATE (comm_j(ng)%jzSend(nMess_j*2))
       ALLOCATE (comm_j(ng)%procRecv(nMess_j*2))
       ALLOCATE (comm_j(ng)%iaRecv(nMess_j*2))
       ALLOCATE (comm_j(ng)%izRecv(nMess_j*2))
       ALLOCATE (comm_j(ng)%jaRecv(nMess_j*2))
       ALLOCATE (comm_j(ng)%jzRecv(nMess_j*2))
    END DO
    
     !X direction
     nsend_i=0
     DO ng=1,ngrids
	i=1
	istep=1
	DO WHILE (i<=nnxp(ng))
	   i=i+istep
	   DO j=2,nnyp(ng)
	      !If there are difference between this point and the last one, we are at border
	       IF (procMap(ng,i,j)/=procMap(ng,i,j-1) .and. procMap(ng,i,j-1)/=0 &
	            .and. procMap(ng,i,j)/=0) THEN
                  p1=procMap(ng,i,j)!;print *,'p1=',p1
	          p2=procMap(ng,i,j-1)!;print *,'p2=',p2
		  IF(myNum==p1) west(ng)=.true.
		  IF(myNum==p2) east(ng)=.true.
	          ini=i-1 !get the first point
	          !There is communication 
		  !p1 -> p2
		  !p2 -> p1	   
	          !looking for until the communication ends. When p1 or p2 differs.
	          DO ii=i+1,nnxp(ng)
	            IF(procMap(ng,ii,j)/=p1 .or. procMap(ng,ii,j-1)/=p2) THEN
		       istep=ii-ini+1
		       EXIT
	            END IF
	            fin=ii+1 !get the last point
	          END DO !ii loop 
	   
	          tagcount=tagcount+1
!GhostZoneLength	   
	          igz=0
		  !makes the comm p1->p2
		  IF(procMap(ng,i,j)==myNum .or. procMap(ng,i,j-1)==myNum) THEN
	             nsend_i=nSend_i+1
	             !WRITE(20+mynum,FMT='("p1->p2: ",7(I2.2,1X))') &
	             !	 p1,p2,ini,fin, &
	             !	 ini-nodei0(p1,ng),fin-nodei0(p1,ng)
	             comm_i(ng)%tag(nsend_i)=tagcount
		     IF(west(ng) .or. east(ng)) igz=igz+GhostZoneLength
	             comm_i(ng)%length(nsend_i)=(fin-ini+1)*igz
	             !
		     comm_i(ng)%procSend(nsend_i)=p1
		     !
	             comm_i(ng)%iaSend(nsend_i)=ini-nodei0(p1,ng)
	             comm_i(ng)%izSend(nsend_i)=fin-nodei0(p1,ng)
!LFR> 	             comm_i(ng)%jaSend(nsend_i)=j  -nodej0(p1,ng)-1+igz	     
!LFR> 	             comm_i(ng)%jzSend(nsend_i)=j  -nodej0(p1,ng)-1+2*igz-1
	             comm_i(ng)%jaSend(nsend_i)=igz+1   
	             comm_i(ng)%jzSend(nsend_i)=2*igz
	             !
		     comm_i(ng)%procRecv(nsend_i)=p2
		     !
	             comm_i(ng)%iaRecv(nsend_i)=ini-nodei0(p2,ng)
	             comm_i(ng)%izRecv(nsend_i)=fin-nodei0(p2,ng)
	             comm_i(ng)%jaRecv(nsend_i)=j  -nodej0(p2,ng)    
	             comm_i(ng)%jzRecv(nsend_i)=j  -nodej0(p2,ng)-1+igz
	          END IF !send i for p1
	   
	          tagcount=tagcount+1
	          igz=0	   
	          !makes the comm p2->p1
	          IF(procMap(ng,i,j)==myNum .or. procMap(ng,i,j-1)==myNum) THEN
	             nSend_i=nSend_i+1	   
		     !WRITE(20+mynum,FMT='("p2->p1: ",7(I2.2,1X))') &
	            !	 p2,p1,ini,fin, &
	            !	 ini-nodei0(p2,ng),fin-nodei0(p2,ng)
	             comm_i(ng)%tag(nSend_i)=tagcount
		     IF(west(ng) .or. east(ng)) igz=igz+GhostZoneLength
	             comm_i(ng)%length(nSend_i)=(fin-ini+1)*igz
	             !
		     comm_i(ng)%procSend(nSend_i)=p2
	             !
		     comm_i(ng)%iaSend(nSend_i)=ini-nodei0(p2,ng)
	             comm_i(ng)%izSend(nSend_i)=fin-nodei0(p2,ng)
	             comm_i(ng)%jaSend(nSend_i)=j   -nodej0(p2,ng)-igz	     
	             comm_i(ng)%jzSend(nSend_i)=j   -nodej0(p2,ng)-1	     
	             !
		     comm_i(ng)%procRecv(nSend_i)=p1
	             !
		     comm_i(ng)%iaRecv(nSend_i)=ini-nodei0(p1,ng)
	             comm_i(ng)%izRecv(nSend_i)=fin-nodei0(p1,ng)
!LFR> 	             comm_i(ng)%jaRecv(nSend_i)=j    -nodej0(p1,ng)-igz+1
!LFR> 	             comm_i(ng)%jzRecv(nSend_i)=j    -nodej0(p1,ng)
	             comm_i(ng)%jaRecv(nSend_i)=1
	             comm_i(ng)%jzRecv(nSend_i)=igz
	          END IF  !send i for p2
	          			      
	        END IF !found Diferences
	   END DO !j loop
        END DO ! i loop
     END DO ! ngrids loop	
	   
     
     nsend_j=0
     DO ng=1,ngrids
	j=1
	jstep=1
	DO WHILE (j<=nnyp(ng))
	   j=j+jstep
	   DO i=2,nnxp(ng)
	      !Se existe diferenca entre esse ponto e o anterior
	      IF (procMap(ng,i,j)/=procMap(ng,i-1,j) .and. procMap(ng,i-1,j)/=0 &
	         .and. procMap(ng,i,j)/=0) THEN
	         p1=procMap(ng,i,j)!;print *,'p1=',p1
	         p2=procMap(ng,i-1,j)!;print *,'p2=',p2
		 IF(myNum==p1) north(ng)=.true.
		 IF(myNum==p2) south(ng)=.true.
	         ini=j-1!;IF(procMap(ng,i,j)==2) PRINT *,'INI=',ini
	         !existe comunicacao entre:
		 !p1 -> p2
		 !p2 -> p1	   
	         !verificando ateh onde p1 fala com p2
	         DO jj=j+1,nnyp(ng)
	           !Se um dos dois eh diferente acabou
	           IF(procMap(ng,i,jj)/=p1 .or. procMap(ng,i-1,jj)/=p2) THEN
		      jstep=jj-ini+1
		      EXIT
	           END IF
	           fin=jj+1
	         END DO 
	   
	         tagcount=tagcount+1
                 igz=0
	         IF(procMap(ng,i,j)==myNum .or. procMap(ng,i-1,j)==myNum) THEN
	            nsend_j=nSend_j+1
	            !WRITE(20+mynum,FMT='("p1->p2: ",7(I2.2,1X))') &
		    !   p1,p2,ini,fin, &
		    !  ini-nodej0(p1,ng),fin-nodej0(p1,ng)
	            comm_j(ng)%tag(nsend_j)=tagcount
		     IF(south(ng) .or. north(ng)) igz=igz+GhostZoneLength
	            comm_j(ng)%length(nsend_j)=(fin-ini+1)*igz
	            comm_j(ng)%procSend(nsend_j)=p1
!LFR> 	            comm_j(ng)%iaSend(nsend_j)=i   -nodei0(p1,ng)-1+igz
!LFR> 	            comm_j(ng)%izSend(nsend_j)=i   -nodei0(p1,ng)-1+2*igz-1
	            comm_j(ng)%iaSend(nsend_j)=igz+1
	            comm_j(ng)%izSend(nsend_j)=2*igz
	            comm_j(ng)%jaSend(nsend_j)=ini-nodej0(p1,ng)	  
	            comm_j(ng)%jzSend(nsend_j)=fin-nodej0(p1,ng)	  
	            comm_j(ng)%procRecv(nsend_j)=p2
   	            comm_j(ng)%iaRecv(nsend_j)=i  -nodei0(p2,ng)
	            comm_j(ng)%izRecv(nsend_j)=i  -nodei0(p2,ng)-1+igz
	            comm_j(ng)%jaRecv(nsend_j)=ini-nodej0(p2,ng)	  
	            comm_j(ng)%jzRecv(nsend_j)=fin-nodej0(p2,ng)
	       END IF

               igz=0
	       tagcount=tagcount+1
	       IF(procMap(ng,i,j)==myNum .or. procMap(ng,i-1,j)==myNum) THEN	   
	          nSend_j=nsend_j+1	   
		  ! WRITE(20+mynum,FMT='("p2->p1: ",7(I2.2,1X))') &
		  !p2,p1,ini,fin, &
		  !ini-nodej0(p2,ng),fin-nodej0(p2,ng)
	         comm_j(ng)%tag(nsend_j)=tagcount
		 IF(south(ng) .or. north(ng)) igz=igz+GhostZoneLength
	         comm_j(ng)%length(nsend_j)=(fin-ini+1)*igz
	         comm_j(ng)%procSend(nsend_j)=p2
	    comm_j(ng)%iaSend(nsend_j)=i  -nodei0(p2,ng)-igz
	    comm_j(ng)%izSend(nsend_j)=i  -nodei0(p2,ng)-1
	         comm_j(ng)%jaSend(nsend_j)=ini-nodej0(p2,ng)	  
	         comm_j(ng)%jzSend(nsend_j)=fin-nodej0(p2,ng)	  
	         comm_j(ng)%procRecv(nsend_j)=p1
!LFR> 	         comm_j(ng)%iaRecv(nsend_j)=i  -nodei0(p1,ng)-igz+1
!LFR> 	         comm_j(ng)%izRecv(nsend_j)=i  -nodei0(p1,ng)
	         comm_j(ng)%iaRecv(nsend_j)=1
	         comm_j(ng)%izRecv(nsend_j)=igz
	         comm_j(ng)%jaRecv(nsend_j)=ini-nodej0(p1,ng)	  
	         comm_j(ng)%jzRecv(nsend_j)=fin-nodej0(p1,ng)
	      END IF  
	       
	  END IF
       END DO
   END DO
   END DO

  END SUBROUTINE PrepareComm

  subroutine CommAlloc(myId,ngrid,m1,m2,m3,ia,iz,ja,jz)
      INTEGER,INTENT(IN) :: myId,ngrid,m1,m2,m3,ia,iz,ja,jz
      
      integer :: i,j,nRec_i,nSnd_i,nRec_j,nSnd_j
      ALLOCATE(newM2(ngrids),newM3(ngrids),newIa(ngrids),newJa(ngrids))
   
     !Counting for allocation
     nRec_i=0
     nSnd_i=0
     nRec_j=0
     nSnd_j=0
     DO i=1,nSend_i
         IF(comm_i(ngrid)%procRecv(i)==myId) nRec_i=nRec_i+1
         IF(comm_i(ngrid)%procSend(i)==myId) nSnd_i=nSnd_i+1
     END DO
     DO j=1,nSend_j
         IF(comm_j(ngrid)%procRecv(j)==myId) nRec_j=nRec_j+1
         IF(comm_j(ngrid)%procSend(j)==myId) nSnd_j=nSnd_j+1
     END DO
 
     !    I Receive
     ALLOCATE(procRecv_i(nRec_i))      
     ALLOCATE(tagRecv_i(nRec_i))       
     ALLOCATE(iaRecv_i(nRec_i))        
     ALLOCATE(izRecv_i(nRec_i))        
     ALLOCATE(jaRecv_i(nRec_i))        
     ALLOCATE(jzRecv_i(nRec_i))        
     ALLOCATE(bufRecvLength_i(nRec_i)) 
     ALLOCATE(bufRecvStart_i(nRec_i))  
     !    I Send		       
     ALLOCATE(procSend_i(nSnd_i))      
     ALLOCATE(tagSend_i(nSnd_i))       
     ALLOCATE(iaSend_i(nSnd_i))        
     ALLOCATE(izSend_i(nSnd_i))        
     ALLOCATE(jaSend_i(nSnd_i))        
     ALLOCATE(jzSend_i(nSnd_i))        
     ALLOCATE(bufSendLength_i(nSnd_i)) 
     ALLOCATE(bufSendStart_i(nSnd_i))  
     !    J Receive		       
     ALLOCATE(procRecv_j(nRec_j))      
     ALLOCATE(tagRecv_j(nRec_j))       
     ALLOCATE(iaRecv_j(nRec_j))        
     ALLOCATE(izRecv_j(nRec_j))        
     ALLOCATE(jaRecv_j(nRec_j))        
     ALLOCATE(jzRecv_j(nRec_j))        
     ALLOCATE(bufRecvLength_j(nRec_j)) 
     ALLOCATE(bufRecvStart_j(nRec_j))  
     !    J Send		       
     ALLOCATE(procSend_j(nSnd_j))      
     ALLOCATE(tagSend_j(nSnd_j))       
     ALLOCATE(iaSend_j(nSnd_j))        
     ALLOCATE(izSend_j(nSnd_j))        
     ALLOCATE(jaSend_j(nSnd_j))        
     ALLOCATE(jzSend_j(nSnd_j))        
     ALLOCATE(bufSendLength_j(nSnd_j)) 
     ALLOCATE(bufSendStart_j(nSnd_j))     
     
     !Local communication and advection buffer with ghostzone
     newM2(ngrid)=iz
     newM3(ngrid)=jz
     newIa(ngrid)=1
     newJa(ngrid)=1

     IF(north(ngrid))  THEN
        newIa(ngrid)=ghostZoneLength+1
     ELSE
        newia(ngrid)=ia 
     END IF
     IF(west(ngrid)) THEN
         newJa(ngrid)=ghostZoneLength+1
     ELSE
        newja(ngrid)=ja
     END IF

     IF(west(ngrid))  newM3(ngrid)=newM3(ngrid)+ghostZoneLength
     IF(east(ngrid))  newM3(ngrid)=newM3(ngrid)+ghostZoneLength
     IF(.not. west(ngrid) .and. .not. east(ngrid)) newM3(ngrid)=m3

     IF(north(ngrid)) newM2(ngrid)=newM2(ngrid)+ghostZoneLength
     IF(south(ngrid)) newM2(ngrid)=newM2(ngrid)+ghostZoneLength
     IF(.not. north(ngrid) .and. .not. south(ngrid)) newM2(ngrid)=m2


      WRITE (70+myId,FMT='(2(I2.2,1X),"N,S,L,W:",4(L1,1X),4(I2.2,1X))') &
		 m2,m3,north(ngrid),south(ngrid),east(ngrid),west(ngrid), &
	   newM2(ngrid),newM3(ngrid),newIa(ngrid),newJa(ngrid)
      CALL flush(70+myId)
     
     
     ALLOCATE(buffcomm(m1,newM2(ngrid),newM3(ngrid)))

  end subroutine CommAlloc

  subroutine CommDeAlloc()
 
     !    I Receive
     DEALLOCATE(procRecv_i)     
     DEALLOCATE(tagRecv_i)	
     DEALLOCATE(iaRecv_i)	
     DEALLOCATE(izRecv_i)	
     DEALLOCATE(jaRecv_i)	
     DEALLOCATE(jzRecv_i)	
     DEALLOCATE(bufRecvLength_i)
     DEALLOCATE(bufRecvStart_i) 
     !    I Send
     DEALLOCATE(procSend_i)	
     DEALLOCATE(tagSend_i)	
     DEALLOCATE(iaSend_i)	
     DEALLOCATE(izSend_i)	
     DEALLOCATE(jaSend_i)	
     DEALLOCATE(jzSend_i)	
     DEALLOCATE(bufSendLength_i)
     DEALLOCATE(bufSendStart_i) 
     !    J Receive		
     DEALLOCATE(procRecv_j)	
     DEALLOCATE(tagRecv_j)	
     DEALLOCATE(iaRecv_j)	
     DEALLOCATE(izRecv_j)	
     DEALLOCATE(jaRecv_j)	
     DEALLOCATE(jzRecv_j)	
     DEALLOCATE(bufRecvLength_j)
     DEALLOCATE(bufRecvStart_j) 
     !    J Send
     DEALLOCATE(procSend_j)	
     DEALLOCATE(tagSend_j)	
     DEALLOCATE(iaSend_j)	
     DEALLOCATE(izSend_j)	
     DEALLOCATE(jaSend_j)	
     DEALLOCATE(jzSend_j)	
     DEALLOCATE(bufSendLength_j)
     DEALLOCATE(bufSendStart_j) 
     DEALLOCATE(buffcomm)
     
  end subroutine CommDeAlloc

   subroutine BuildComm(myId,ngrid,m1,nRec_i,nSnd_i,nRec_j,nSnd_j, &
                        bufSendTotalLength_i,bufSendTotalLength_j, &
			bufRecvTotalLength_i,bufRecvTotalLength_j)
	   
   integer, intent(in)  :: myId,ngrid,m1
   integer, intent(out) :: nRec_i,nSnd_i,nRec_j,nSnd_j
   integer, intent(out) :: bufSendTotalLength_i,bufSendTotalLength_j
   integer, intent(out) :: bufRecvTotalLength_i,bufRecvTotalLength_j

   integer :: i,j
         	   
  nRec_i=0
  nSnd_i=0
  nRec_j=0
  nSnd_j=0
  bufSendTotalLength_i=0
  bufSendTotalLength_j=0
  bufRecvTotalLength_i=0
  bufREcvTotalLength_j=0
  DO i=1,nSend_i
      IF(comm_i(ngrid)%procRecv(i)==myId) THEN
         nRec_i=nRec_i+1
         bufRecvStart_i(nRec_i)=bufRecvTotalLength_i+1
         procRecv_i(nRec_i)=comm_i(ngrid)%procSend(i)
         tagRecv_i(nRec_i)=comm_i(ngrid)%tag(i)
         iaRecv_i(nRec_i)=comm_i(ngrid)%iaRecv(i)
         izRecv_i(nRec_i)=comm_i(ngrid)%izRecv(i)
         jaRecv_i(nRec_i)=comm_i(ngrid)%jaRecv(i)
         jzRecv_i(nRec_i)=comm_i(ngrid)%jzRecv(i)
         bufRecvLength_i(nRec_i)=comm_i(ngrid)%length(i)*m1
	 bufRecvTotalLength_i=bufRecvTotalLength_i+bufRecvLength_i(nRec_i)
	 
	 !
      END IF
      IF(comm_i(ngrid)%procSend(i)==myId) THEN
         nSnd_i=nSnd_i+1
         bufSendStart_i(nSnd_i)=bufSendTotalLength_i+1
         procSend_i(nSnd_i)=comm_i(ngrid)%procRecv(i)
         tagSend_i(nSnd_i)=comm_i(ngrid)%tag(i)
         iaSend_i(nSnd_i)=comm_i(ngrid)%iaSend(i)
         izSend_i(nSnd_i)=comm_i(ngrid)%izSend(i)
         jaSend_i(nSnd_i)=comm_i(ngrid)%jaSend(i)
         jzSend_i(nSnd_i)=comm_i(ngrid)%jzSend(i)
         bufSendLength_i(nSnd_i)=comm_i(ngrid)%length(i)*m1
	 bufSendTotalLength_i=bufSEndTotalLength_i+bufSendLength_i(nSnd_i)
      END IF
  END DO
  DO j=1,nSend_j
      IF(comm_j(ngrid)%procRecv(j)==myId) THEN
         nRec_j=nRec_j+1
         bufRecvStart_j(nRec_j)=bufRecvTotalLength_j+1
         procRecv_j(nRec_j)=comm_j(ngrid)%procSend(j)
         tagRecv_j(nRec_j)=comm_j(ngrid)%tag(j)
         iaRecv_j(nRec_j)=comm_j(ngrid)%iaRecv(j)
         izRecv_j(nRec_j)=comm_j(ngrid)%izRecv(j)
         jaRecv_j(nRec_j)=comm_j(ngrid)%jaRecv(j)
         jzRecv_j(nRec_j)=comm_j(ngrid)%jzRecv(j)
         bufRecvLength_j(nRec_j)=comm_j(ngrid)%length(j)*m1
	 bufRecvTotalLength_j=bufrEcvTotalLength_j+bufRecvLength_j(nRec_j)
      END IF
      IF(comm_j(ngrid)%procSend(j)==myId) THEN
         nSnd_j=nSnd_j+1
         bufSendStart_j(nSnd_j)=bufSendTotalLength_j+1
         procSend_j(nSnd_j)=comm_j(ngrid)%procRecv(j)
         tagSend_j(nSnd_j)=comm_j(ngrid)%tag(j)
         iaSend_j(nSnd_j)=comm_j(ngrid)%iaSend(j)
         izSend_j(nSnd_j)=comm_j(ngrid)%izSend(j)
         jaSend_j(nSnd_j)=comm_j(ngrid)%jaSend(j)
         jzSend_j(nSnd_j)=comm_j(ngrid)%jzSend(j)
         bufSendLength_j(nSnd_j)=comm_j(ngrid)%length(j)*m1
	 bufSendTotalLength_j=bufSendTotalLength_j+bufSendLength_j(nSnd_j)
      END IF
   END DO
 
   end subroutine BuildComm

   subroutine DumpBorders(prefix,dir,m1, m2, m3, mynum, &
	nRcv, procRecv, tagRecv, iaRecv, izRecv, jaRecv, jzRecv, &
	bufRecvStart, bufRecvLength, bufRecvTotalLength, &
	nSnd, procSend, tagSend, iaSend, izSend, jaSend, jzSend, &
	bufSendStart, bufSendLength, bufSendTotalLength)
     
     CHARACTER(len=*),intent(in) :: prefix
     CHARACTER(len=1),intent(in) :: dir
     integer, intent(in) :: m1
     integer, intent(in) :: m2
     integer, intent(in) :: m3,mynum
     integer, intent(in) :: nRcv
     integer, intent(in) :: procRecv(nRcv)
     integer, intent(in) :: tagRecv(nRcv)
     integer, intent(in) :: iaRecv(nRcv)
     integer, intent(in) :: izRecv(nRcv)
     integer, intent(in) :: jaRecv(nRcv)
     integer, intent(in) :: jzRecv(nRcv)
     integer, intent(in) :: bufRecvStart(nRcv)
     integer, intent(in) :: bufRecvLength(nRcv)
     integer, intent(in) :: bufRecvTotalLength
     integer, intent(in) :: nSnd
     integer, intent(in) :: procSend(nSnd)
     integer, intent(in) :: tagSend(nSnd)
     integer, intent(in) :: iaSend(nSnd)
     integer, intent(in) :: izSend(nSnd)
     integer, intent(in) :: jaSend(nSnd)
     integer, intent(in) :: jzSend(nSnd)
     integer, intent(in) :: bufSendStart(nSnd)
     integer, intent(in) :: bufSendLength(nSnd)
     integer, intent(in) :: bufSendTotalLength
 
     INTEGER :: i,unitout
     CHARACTER(LEN=2) :: c_mynum
     
     unitout=40+mynum
     WRITE(c_mynum,FMT='(I2.2)') mynum
     
     OPEN (UNIT=unitout,file=prefix//dir//'-'//c_mynum//'.dat')

     WRITE (unitout,FMT='("=== Receive === ",I3.3," ===",4(I3.3,1X),I8') mynum,m1,m2,m3,nRcv,bufRecvTotalLength 
     CALL flush(unitout)
        
     WRITE (unitout,FMT='(A3,1X,A3,1X,A8,1X,4(A3,1X),A8,1X,A8)') &
	                  'Nrc','Prc',' Tag ','ia','iz','ja','jz','start','lenght';CALL flush(unitout)
     DO i=1,nRcv
        WRITE (unitout,FMT='(I3.3,1X,I3.3,1X,I8.8,1X,4(I3.3,1X),I8,1X,I8)') &
	      i,procRecv(i),tagRecv(i),iaRecv(i),izRecv(i),jaRecv(i), &
	      jzRecv(i),bufRecvStart(i),bufRecvLength(i);CALL flush(unitout)
     END DO

     WRITE (unitout,FMT='("===   Send  === ",I3.3," ===",4(I3.3,1X),I8') mynum,m1,m2,m3,nSnd,bufSendTotalLength    
     WRITE (unitout,FMT='(A3,1X,A3,1X,A8,1X,4(A3,1X),A8,1X,A8)') &
	                  'Nsn','Prc',' Tag ','ia','iz','ja','jz','start','lenght';CALL flush(unitout)
     DO i=1,nSnd
        WRITE (unitout,FMT='(I3.3,1X,I3.3,1X,I8.8,1X,4(I3.3,1X),I8,1X,I8)') &
	      i,procSend(i),tagSend(i),iaSend(i),izSend(i),jaSend(i), &
	      jzSend(i),bufSendStart(i),bufSendLength(i);CALL flush(unitout)
     END DO
     CALL flush(unitout)

     CLOSE(unitout)
     
  END SUBROUTINE DumpBorders
  
  SUBROUTINE InitialFieldsUpdate(m1,m2,m3,Nm2,Nm3,ng,mynum, &
	nRec, procRecv, tagRecv, iaRecv, izRecv, jaRecv, jzRecv, &
	bufRecvStart, bufRecvLength, bufRecvTotalLength, &
	nSnd, procSend, tagSend, iaSend, izSend, jaSend, jzSend, &
	bufSendStart, bufSendLength, bufSendTotalLength)
  
     integer, intent(in) :: m1
     integer, intent(in) :: m2,nm2
     integer, intent(in) :: m3,nm3,ng,mynum
     integer, intent(in) :: nRec
     integer, intent(in) :: procRecv(nRec)
     integer, intent(in) :: tagRecv(nRec)
     integer, intent(in) :: iaRecv(nRec)
     integer, intent(in) :: izRecv(nRec)
     integer, intent(in) :: jaRecv(nRec)
     integer, intent(in) :: jzRecv(nRec)
     integer, intent(in) :: bufRecvStart(nRec)
     integer, intent(in) :: bufRecvLength(nRec)
     integer, intent(in) :: bufRecvTotalLength
     integer, intent(in) :: nSnd
     integer, intent(in) :: procSend(nSnd)
     integer, intent(in) :: tagSend(nSnd)
     integer, intent(in) :: iaSend(nSnd)
     integer, intent(in) :: izSend(nSnd)
     integer, intent(in) :: jaSend(nSnd)
     integer, intent(in) :: jzSend(nSnd)
     integer, intent(in) :: bufSendStart(nSnd)
     integer, intent(in) :: bufSendLength(nSnd)
     integer, intent(in) :: bufSendTotalLength

     REAL,DIMENSION (m1,nm2,nm3) :: l_dxtW,l_dytW
     INTEGER :: i,j,k
     
    IF(bufSendTotalLength==0 .or. bufRecvTotalLength==0) RETURN
    
     DO i=1,nm2
     	DO j=1,nm3
     	   DO k=1,m1
     	      l_dxtW(k,i,j)=advmntNew_g(ng)%dxtW(i,j)
     	      l_dytW(k,i,j)=advmntNew_g(ng)%dytW(i,j)
     	   END DO
     	END DO
     END DO


    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%u3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,1; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%v3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,2; CALL flush(6)
     
    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%w3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,3; CALL flush(6)
    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%vc3d_jn, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,4; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%vc3d_out, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,5; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%vc3d_x, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,6; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%vc3d_y, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,7; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%dd0_3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,8; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%dd0_3du, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,9; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%dd0_3dv, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,10; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%dd0_3dw, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,11; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%den0_3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,12; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%den1_3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,13; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%den2_3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,14; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%den3_3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,15; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,l_dxtW, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,16; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,l_dytW, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,17; CALL flush(6)

      DO i=1,nm2
     	DO j=1,nm3
     	      advmntNew_g(ng)%dxtW(i,j)=l_dxtW(1,i,j)
     	      advmntNew_g(ng)%dytW(i,j)=l_dytW(1,i,j)
     	END DO
     END DO

  END SUBROUTINE InitialFieldsUpdate

  SUBROUTINE CopyInitFields(ngrid,m1,m2,m3)
     INTEGER, INTENT(IN) :: ngrid,m1,m2,m3
     INTEGER :: iBegin,iEnd,jBegin,jEnd
     
     iBegin=newIa(ngrid)-1
     iEnd=iBegin+m2-1
     jBegin=newJa(ngrid)-1
     jEnd=jBegin+m3-1
     
!     IF(bigDump==1) 
   WRITE (*,fmt='(A,5(I3.3,1X))') 'Copy limits: ',mynum,iBegin,iEnd,jBegin,jEnd
     
   advmntNew_g(ngrid)%u3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%     u3d(1:m1,1:m2,1:m3) 
   advmntNew_g(ngrid)%v3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%     v3d(1:m1,1:m2,1:m3) 
   advmntNew_g(ngrid)%w3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%     w3d(1:m1,1:m2,1:m3) 

   !advmntNew_g(ngrid)%vc3d_jn(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% vc3d_jn(1:m1,1:m2,1:m3)
   !advmntNew_g(ngrid)%vc3d_out(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%vc3d_out(1:m1,1:m2,1:m3)
   advmntNew_g(ngrid)%vc3d_x(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%  vc3d_x(1:m1,1:m2,1:m3) 
   advmntNew_g(ngrid)%vc3d_y(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%  vc3d_y(1:m1,1:m2,1:m3) 

   advmntNew_g(ngrid)%dd0_3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%  dd0_3d(1:m1,1:m2,1:m3) 
   advmntNew_g(ngrid)%dd0_3du(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% dd0_3du(1:m1,1:m2,1:m3)
   advmntNew_g(ngrid)%dd0_3dv(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% dd0_3dv(1:m1,1:m2,1:m3)
   advmntNew_g(ngrid)%dd0_3dw(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% dd0_3dw(1:m1,1:m2,1:m3)

   advmntNew_g(ngrid)%den0_3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% den0_3d(1:m1,1:m2,1:m3)
   advmntNew_g(ngrid)%den1_3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% den1_3d(1:m1,1:m2,1:m3)
   advmntNew_g(ngrid)%den2_3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% den2_3d(1:m1,1:m2,1:m3)
   advmntNew_g(ngrid)%den3_3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% den3_3d(1:m1,1:m2,1:m3)

   advmntNew_g(ngrid)%dxtW(iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%    dxtW(1:m2,1:m3) 
   advmntNew_g(ngrid)%dytW(iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%    dytW(1:m2,1:m3) 
   advmntNew_g(ngrid)%dztW(1:m1)=advmnt_g(ngrid)%    dztW(1:m1)   

  END SUBROUTINE CopyInitFields

 
SUBROUTINE Advec3d_X(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn,mynum)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !

  use node_mod, only : nodeibcon,nodeia,nodeiz,nodeja,nodejz, &
                       ipaths,iget_paths,nodei0,nodej0,nmachs
       
  use mem_grid, only:  ngrids, nnxp,nnyp

  IMPLICIT none  
 
  INCLUDE "i8.h"
 
  INTEGER,INTENT(IN)                        :: m1,m2,m3, ia,iz,ja,jz,mynum
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)    :: u
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)    :: den0
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)    :: den1
  REAL   ,INTENT(in)                        :: dt
  REAL   ,INTENT(in),DIMENSION(m2,m3)       :: dxx
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)    :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)   :: qn

  REAL,DIMENSION(m1,m2,m3)    :: flux
  REAL,DIMENSION(m1,m2,m3)    :: vcmax
  REAL,DIMENSION(m1,m2,m3)    :: vcmin
  LOGICAL,DIMENSION(m1,m2,m3) :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: idime
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  
  INTEGER :: ii,ji,ii0,ji0,ie,je,ie0,je0,ipos,iia,iiz,nvar
  INTEGER,DIMENSION(nmachs*nmachs*ngrids*2*nvariables) :: isend_req,irecv_req
  
  
  integer :: nf
  
  !nvariables=
  !             1- flux
  REAL,DIMENSION(nvariables,m1) :: varSend
  REAL,DIMENSION(nvariables,m1) :: varRecv  
  
  nf=80+mynum
 idime = m2 ! x dir

  qn=q0

! WRITE (nf,FMT='("Idime= ",I3.3)') idime
! WRITE (nf,FMT='("m2=",I3.3," ia=",I3.3," iz=",I3.3," Map(m2)=",I3.3," Map(ia)=",I3.3," Map(iz)=",I3.3)') &
!       m2,ia,iz,m2+nodei0(mynum,ngrid),ia+nodei0(mynum,ngrid),iz+nodei0(mynum,ngrid);CALL flush(nf)
 
  imxmn=.false.
  ! imxmn(idime+1)=.false.

 DO j=ja,jz

    DO k=2,m1-1 
       ! Update mixing ratios and limit Fluxes going UP where u>0
       !  First assume upstream flux at edge of domain
       IF(u(k,1,j)>=zr0) flux(k,1,j)= q0(k,1,j)*u(k,1,j)*dt*dd0(k,1,j)
    END DO
  END DO 
 
 DO j=ja,jz
 
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  i=2,idime-1 ! ia,iz-1 or 1,iz-1

     DO k=2,m1-1 
				  
        imxmn(k,i,j)=q0(k,i,j)>=max(q0(k,i-1,j),q0(k,i+1,j)) .or. & !=true if local
              q0(k,i,j)<=min(q0(k,i-1,j),q0(k,i+1,j))       !       extrema
        ck1= q0(k,i,j)
        ck2= q0(k,i,j)
        if(u(k,i,j  )< zr0) ck1= q0(k,i+1,j)
        if(u(k,i-1,j)>=zr0) ck2= q0(k,i-1,j)
        vcmax(k,i,j)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        vcmin(k,i,j)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        !  VCMAX and VCMIN are the absolute physical limits to the
        !     mixing ratio at t+dt. If these limits are ever violated,
        !     non-monotonic (oscillatory) behavior in solution results
        !

        IF(u(k,i,j)<zr0) CYCLE
     
        IF(u(k,i-1,j)<zr0) THEN
           flux(k,i,j)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell     
        ELSE                              !      use upstream
           x1= dt*u(k,i,j)/dxx(i,j)               ! Courant number
           x1n= (1.-x1)*(q0(k,i+1,j)-q0(k,i-1,j))/4.
           !
           ! First, estimate mixing ratio in outflowing fluid (Cf)
           cf= q0(k,i,j) + x1n                                       !Eq-4a
           !
           !   Check to see if there is a peak (min) upwind and/or
           !    downwind of cell face
           IF(imxmn(k,i-1,j)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(k,i+1,j)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
           !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
           !
           !   Limit Cf to be between mixing ratio on either side of edge
           !      where flux is being calculated
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i+1,j))  ), MAX(q0(k,i,j),q0(k,i+1,j)) )
           !- for debug purposes only
	   if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
           !                       implemented this subroutine constant mixing
           !                       ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	   !
           !   Calculate mixing ratio at new time, but limit to physically
           !    reasonable values
           qn(k,i,j)= MAX(vcmin(k,i,j),MIN(vcmax(k,i,j),          &   !eq-3&8
                    (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(k,i-1,j)/dxx(i,j))/den1(k,i,j) ))

	   ! for debug purposes only
           if(debug==on)  QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-X1*CF1*DD0(k,i,j)+FLUX(k,I-1,j)/DXX(I,J))/DEN1(k,i,j)
           !            This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
	   !                       implemented this subroutine constant mixing
           !                       ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
           !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	   !       if(qn(i) < 99.999) pause

	   !
           !   Re-calculate OUTFLOWING flux before moving on to next cell
           !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
           !    is encountered.
           flux(k,i,j)= dxx(i,j)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(k,i-1,j)
        END IF                                                  !Eq-9a
     END DO
     
  END DO
END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !

 DO j=ja,jz

    DO k=2,m1-1       
       IF(u(k,idime-1,j)<zr0) flux(k,idime-1,j)= &
    	       q0(k,idime,j)*u(k,idime-1,j)*dt*dd0(k,idime-1,j)
    END DO
  END DO 


  DO j=ja,jz
   
   !LFR-> from iz to ia
   DO i=idime-1,2,-1 !iz,ia,-1
     ! print*,'part 2'
     DO k=2,m1-1 
        
        IF(u(k,i-1,j)>=zr0) THEN           ! Inflow-only cell
        
	   IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(k,i,j),   MIN(   vcmax(k,i,j),&
                  (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j) + flux(k,i-1,j)/dxx(i,j))/den1(k,i,j) ))
           ELSE
              x1=  dt*ABS(u(k,i-1,j))/dxx(i,j)     ! Courant number
              x1n= (1.-x1)*(q0(k,i-1,j)-q0(k,i+1,j))/4.
              cf= q0(k,i,j) + x1n                                       !Eq-4b
              IF(imxmn(k,i+1,j)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
              IF(imxmn(k,i-1,j)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
              cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i-1,j)) ), MAX(q0(k,i,j),q0(k,i-1,j)) )
              !- for debug purposes only
	      if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
              !			   ratios should be maintained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	      IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
              qn(k,i,j)= MAX(  vcmin(k,i,j),  MIN(   vcmax(k,i,j), 	  &   !Eq-3&8
                    (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j)-x1*cf1*dd0(k,i-1,j))/den1(k,i,j) ))
              !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	      !- for debug purposes only
              if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(k,I,j)/DXX(I,J)-X1*CF1*DD0(k,I-1,j))/DEN1(k,i,j)
              !		This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
	      !			   ratios should be maintiained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
              !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
              !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	      !      if(qn(i) < 99.999) pause

	      flux(k,i-1,j)=dxx(i,j)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(k,i,j)!Eq-9b
        END IF                                                  
     END DO

  END DO

  !END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO !- big loop y-z


END SUBROUTINE Advec3d_X

!---------------------------------------------------------------------------

SUBROUTINE Advec3d_Y(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !
  use node_mod, only : nodeibcon,nodeia,nodeiz,nodeja,nodejz, &
                       ipaths,iget_paths,nodei0,nodej0,nmachs
       
  use mem_grid, only:  ngrids, nnxp,nnyp


  IMPLICIT none  
 
  INCLUDE "i8.h"
  
  INTEGER,INTENT(IN)                   :: m1,m2,m3, ia,iz,ja,jz
  INTEGER                  :: idime
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)   :: u
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)     :: den0
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)     :: den1
  REAL   ,INTENT(in)                   :: dt
  REAL   ,INTENT(in),DIMENSION(m2,m3)     :: dxx
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)   :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)    :: qn

  REAL,DIMENSION(m1,m2,m3)                 :: flux
  REAL,DIMENSION(m1,m2,m3)                   :: vcmax
  REAL,DIMENSION(m1,m2,m3)                   :: vcmin
  LOGICAL,DIMENSION(m1,m2,m3)            :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  !
  INTEGER :: ii,ji,ii0,ji0,ie,je,ie0,je0,ipos,iia,iiz,nvar
  INTEGER,DIMENSION(nmachs*nmachs*ngrids*2*nvariables) :: isend_req,irecv_req
  
  
  integer :: nf
  
  !nvariables=
  !             1- q
  !             2- flux
  REAL,DIMENSION(nvariables,m1) :: varSend
  REAL,DIMENSION(nvariables,m1) :: varRecv  
  
  nf=80+mynum

  idime = m3 ! y dir

  qn= q0
  imxmn=.false.
  ! imxmn(idime+1)=.false.

 DO i=ia,iz

    ! Update mixing ratios and limit Fluxes going UP where u>0
    !  First assume upstream flux at edge of domain
    DO k=2,m1-1   
       IF(u(k,i,1)>=zr0) flux(k,i,1)= q0(k,i,1)*u(k,i,1)*dt*dd0(k,i,1)
    END DO
      
  END DO 

 !- big loop y-z
 DO i=ia,iz
   
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  j=2,idime-1 ! ja,jz
     DO k=2,m1-1   
        imxmn(k,i,j)=q0(k,i,j)>=max(q0(k,i,j-1),q0(k,i,j+1)) .or. & !=true if local
                     q0(k,i,j)<=min(q0(k,i,j-1),q0(k,i,j+1))	    !	    extrema
        ck1= q0(k,i,j)
        ck2= q0(k,i,j)
        if(u(k,i,j  )< zr0) ck1= q0(k,i,j+1)
        if(u(k,i,j-1)>=zr0) ck2= q0(k,i,j-1)
        vcmax(k,i,j)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        vcmin(k,i,j)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        !  VCMAX and VCMIN are the absolute physical limits to the
        !	mixing ratio at t+dt. If these limits are ever violated,
        !	non-monotonic (oscillatory) behavior in solution results
        !print*,'part 1'
        IF(u(k,i  ,j)<zr0) CYCLE
     
        IF(u(k,i,j-1)<zr0) THEN
           flux(k,i,j)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell
     
        ELSE                              !      use upstream
           x1= dt*u(k,i,j)/dxx(i,j)               ! Courant number
           x1n= (1.-x1)*(q0(k,i,j+1)-q0(k,i,j-1))/4.
           !
           ! First, estimate mixing ratio in outflowing fluid (Cf)
           cf= q0(k,i,j) + x1n                                       !Eq-4a
           !
           !   Check to see if there is a peak (min) upwind and/or
           !    downwind of cell face
           IF(imxmn(k,i,j-1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(k,i,j+1)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
           !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
           !
           !   Limit Cf to be between mixing ratio on either side of edge
           !      where flux is being calculated
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i,j+1))  ), MAX(q0(k,i,j),q0(k,i,j+1)) )
           !- for debug purposes only
	   if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
           !                       implemented this subroutine constant mixing
           !                       ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	   !
           !   Calculate mixing ratio at new time, but limit to physically
           !    reasonable values
           qn(k,i,j)= MAX(  vcmin(k,i,j),   MIN(   vcmax(k,i,j),          &   !eq-3&8
                 (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(k,i,j-1)/dxx(i,j))/den1(k,i,j) ))

	   ! for debug purposes only
           if(debug==on)  QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-X1*CF1*DD0(k,i,j)+FLUX(k,i,J-1)/DXX(I,J))/DEN1(k,i,j)
           !            This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
	   !                       implemented this subroutine constant mixing
           !                       ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
           !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	   !       if(qn(i) < 99.999) pause

	   !
           !   Re-calculate OUTFLOWING flux before moving on to next cell
           !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
           !    is encountered.
           flux(k,i,j)= dxx(i,j)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(k,i,j-1)
        END IF                                                  !Eq-9a
     END DO
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !
  END DO
  
 DO i=ia,iz

    DO k=2,m1-1   
       IF(u(k,i,idime-1)<zr0) flux(k,i,idime-1)=q0(k,i,idime)*u(k,i,idime-1)*dt*dd0(k,i,idime-1)
    END DO
  END DO 

  DO i=ia,iz
  
  DO j=idime-1,2,-1 !jz,ja,-1
     DO k=2,m1-1     
        IF(u(k,i,j-1)>=zr0) THEN           ! Inflow-only cell        
	   IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(k,i,j),   MIN(   vcmax(k,i,j),&
                  (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j) + flux(k,i,j-1)/dxx(i,j))/den1(k,i,j) ))
           ELSE
              x1=  dt*ABS(u(k,i,j-1))/dxx(i,j)     ! Courant number
              x1n= (1.-x1)*(q0(k,i,j-1)-q0(k,i,j+1))/4.
              cf= q0(k,i,j) + x1n                                       !Eq-4b
              IF(imxmn(k,i,j+1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
              IF(imxmn(k,i,j-1)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
              cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i,j-1)) ), MAX(q0(k,i,j),q0(k,i,j-1)) )
              !- for debug purposes only
	      if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
              !			   ratios should be maintained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	      IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
              qn(k,i,j)= MAX(  vcmin(k,i,j),  MIN(   vcmax(k,i,j), 	  &   !Eq-3&8
                     (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j)-x1*cf1*dd0(k,i,j-1))/den1(k,i,j) ))
              !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	      !- for debug purposes only
              if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(k,i,J)/DXX(I,J)-X1*CF1*DD0(k,I,j-1))/DEN1(k,i,j)
              !		This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
	      !			   ratios should be maintiained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
              !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
              !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	      !      if(qn(i) < 99.999) pause

	      flux(k,i,j-1)=dxx(i,j)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(k,i,j)!Eq-9b
        END IF                                                  
     END DO
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO !- big loop x-z


END SUBROUTINE Advec3d_Y

!---------------------------------------------------------------------------
SUBROUTINE Advec3d_Z(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !
  IMPLICIT none
  
  INTEGER,INTENT(IN)                   :: m1,m2,m3, ia,iz,ja,jz
  INTEGER                  :: idime
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)   :: u
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)     :: den0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)     :: den1
  REAL   ,INTENT(in)                   :: dt
  REAL   ,INTENT(in),DIMENSION(m1)     :: dxx
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)   :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)    :: qn

  REAL,DIMENSION(m1)                 :: flux
  REAL,DIMENSION(m1)                   :: vcmax
  REAL,DIMENSION(m1)                   :: vcmin
  LOGICAL,DIMENSION(m1)            :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  !

 idime = m1 ! z dir

  qn = q0

 !- big loop y-x
 DO j=ja,jz ; DO i=ia,iz
 
  
  imxmn=.false.
  ! imxmn(idime+1)=.false.
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  k=2,idime-1 !
  
     imxmn(k)=q0(k,i,j)>=max(q0(k-1,i,j),q0(k+1,i,j)) .or. & !=true if local
              q0(k,i,j)<=min(q0(k-1,i,j),q0(k+1,i,j))	    !	    extrema
     ck1= q0(k,i,j)
     ck2= q0(k,i,j)
     if(u(k  ,i,j)< zr0) ck1= q0(k+1,i,j)
     if(u(k-1,i,j)>=zr0) ck2= q0(k-1,i,j)
     vcmax(k)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
     vcmin(k)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
  ENDDO
  !  VCMAX and VCMIN are the absolute physical limits to the
  !     mixing ratio at t+dt. If these limits are ever violated,
  !     non-monotonic (oscillatory) behavior in solution results
  !
  ! Update mixing ratios and limit Fluxes going UP where u>0
  !  First assume upstream flux at edge of domain
  IF(u(1,i,j)>=zr0) flux(1)= q0(1,i,j)*u(1,i,j)*dt*dd0(1,i,j)
  
  DO k=2,idime-1 ! ja,jz
  
  !print*,'part 1'
     IF(u(k,i  ,j)<zr0) CYCLE
     
     IF(u(k-1,i,j)<zr0) THEN
        flux(k)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell
     
     ELSE                              !      use upstream
        x1= dt*u(k,i,j)/dxx(k)               ! Courant number
        x1n= (1.-x1)*(q0(k+1,i,j)-q0(k-1,i,j))/4.
        !
        ! First, estimate mixing ratio in outflowing fluid (Cf)
        cf= q0(k,i,j) + x1n                                       !Eq-4a
        !
        !   Check to see if there is a peak (min) upwind and/or
        !    downwind of cell face
        IF(imxmn(k-1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
        IF(imxmn(k+1)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
        !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
        !
        !   Limit Cf to be between mixing ratio on either side of edge
        !      where flux is being calculated
        cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k+1,i,j))  ), MAX(q0(k,i,j),q0(k+1,i,j)) )
        !- for debug purposes only
	if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
        !                       implemented this subroutine constant mixing
        !                       ratios should be maintained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	!
        !   Calculate mixing ratio at new time, but limit to physically
        !    reasonable values
        qn(k,i,j)= MAX(  vcmin(k),   MIN(   vcmax(k),          &   !eq-3&8
        (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(k-1)/dxx(k))/den1(k,i,j) ))

	! for debug purposes only
        if(debug==on)  QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-X1*CF1*DD0(k,i,j)+FLUX(k-1)/dxx(k))/DEN1(k,i,j)
        !            This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
	!                       implemented this subroutine constant mixing
        !                       ratios should be maintiained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
       !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	!       if(qn(i) < 99.999) pause

	!
        !   Re-calculate OUTFLOWING flux before moving on to next cell
        !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
        !    is encountered.
        flux(k)= dxx(k)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(k-1)
     END IF                                                  !Eq-9a
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !
  
  
  IF(u(idime-1,i,j)<zr0) flux(idime-1)=q0(idime,i,j)*u(idime-1,k,j)*dt*dd0(idime-1,i,j)
 
  DO k=idime-1,2,-1 !jz,ja,-1
     
       IF(u(k-1,i,j)>=zr0) THEN           ! Inflow-only cell
        
	IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(k),   MIN(   vcmax(k),&
           (q0(k,i,j)*den0(k,i,j)-flux(k)/dxx(k) + flux(k-1)/dxx(k))/den1(k,i,j) ))
        ELSE
           x1=  dt*ABS(u(k-1,i,j))/dxx(k)     ! Courant number
           x1n= (1.-x1)*(q0(k-1,i,j)-q0(k+1,i,j))/4.
           cf= q0(k,i,j) + x1n                                       !Eq-4b
           IF(imxmn(k+1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(k-1)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k-1,i,j)) ), MAX(q0(k,i,j),q0(k-1,i,j)) )
           !- for debug purposes only
	   if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
           !			   ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	   IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
           qn(k,i,j)= MAX(  vcmin(k),  MIN(   vcmax(k), 	  &   !Eq-3&8
               (q0(k,i,j)*den0(k,i,j)-flux(k)/dxx(k)-x1*cf1*dd0(k-1,i,j))/den1(k,i,j) ))
           !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	   !- for debug purposes only
           if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(k)/dxx(k)-X1*CF1*DD0(k-1,I,j))/DEN1(k,i,j)
           !		This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
	   !			   ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
         !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
         !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
         !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	 !      if(qn(i) < 99.999) pause

	   flux(k-1)=dxx(k)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(k)!Eq-9b
     END IF                                                  
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO;ENDDO !- big loop y-x


END SUBROUTINE Advec3d_Z

     
END MODULE monotonic_adv
